text:
# Observational constraints for Lithium depletion before the RGB.Based on observations collected at ESO, La Silla, Chile, and at the Observatoire de Haute Provence, France, operated by the Centre National de la Recherche Scientifique (CNRS)

## 1 Introduction

Despite the important advances made in the past decade in the study of the stellar lithium behavior, a large number of questions are not yet answered. We do not completely understand the processes controlling lithium production, nor how and when lithium is depleted. Related to these questions are the physical bases of the mixing mechanisms in stellar interiors. Following the initial study by Bonsack (1959) different works have attempted to establish the behavior of lithium abundance along the giant branch (e.g.: Alschuler 1975; Wallerstein 1966; Brown et al. 1989) and its link to rotation (Wallerstein et al. 1994; De Medeiros et al. 2000). These works have shown a steady decline in lithium abundances from spectral types F5III to F8III, a wide spread around G0III and a gradual decline with temperature for stars redward of such spectral type. Wallerstein et al. (1994) have found that giant stars with $`v\mathrm{sin}i`$$`>`$ 50 $`\mathrm{km}\mathrm{s}^1`$ , located in the (B$``$V) color interval from 0.40 to 0.70, present lithium abundance close to the presumed primordial value, whereas slower rotators show reduced lithium abundances in spite of their earlier spectral types. De Medeiros et al. (2000) have found a trend of discontinuity in the distribution of the lithium abundances around the spectral type G0III, paralleling the sudden decline observed in rotational velocity. The origin of such a discontinuity is not yet well established but it seems to strongly depend on stellar mass (De Medeiros and Mayor 1990). Nevertheless, the rotation-lithium connection in giant stars appears to be a more complicated problem. In spite of the fact that high lithium content is associated with fast rotation, slow rotators present a large spread in the values of lithium abundance (De Medeiros et al., 2000). In addition, this connection seems to depend on the stellar mass, metallicity and age. Another important question concerns the level of dilution of lithium along the giant branch. While standard theory predicts a factor of dilution of about 40 to 60 for 1 $`M_{}`$and 2 $`M_{}`$respectively (Iben, 1965a, 1965b), the observations for these stellar masses show that the factor of dilution of lithium is far higher, reaching values as large as 400 to 1000. A more solid discussion about this factor of dilution requires more measurements of lithium abundance for stars located near and at the Hertzsprung gap, in particular for stars blueward of G0III, namely stars at the blue side of the gap. The distribution of lithium abundances before stars evolve along the giant branch may help us to understand the nature of the apparent discontinuity in lithium content around G0III and to establish on a more solid basis the dilution factor of stars evolving along the giant branch.

In this work we present new measurements of lithium abundance for giant stars in the spectral range F5III to G5III, typically stars located near and along the Hertzsprung gap. By combining these data with precise rotational velocities we analyse the rotation-lithium connection and the nature of the discontinuities in the distribution of these two stellar parameters.

## 2 Observations

For this study we have selected 54 luminosity class III giants with spectral types ranging from F5III to G5III. All these stars have been previously observed with the CORAVEL spectrograph (Baranne et al., 1979) for rotational velocity measurements and binarity signatures (De Medeiros and Mayor, 1999).

### 2.1 Rotational velocity

For most of the stars of the present sample, rotational velocity measurements were taken from The Catalog of Rotational and Radial Velocities for Evolved Stars by De Medeiros and Mayor (1999). In this work, the authors present precise rotational velocity, $`v\mathrm{sin}i`$ , for evolved stars, obtained on the basis of observations acquired with the CORAVEL spectrometers (Baranne et al. 1979). For giant stars of luminosity class III in particular, the measured rotational velocities have an uncertainty of about 1.0 $`\mathrm{km}\mathrm{s}^1`$ for stars with $`v\mathrm{sin}i`$ lower than about 30.0 $`\mathrm{km}\mathrm{s}^1`$ . For faster rotators, the estimations indicate an uncertainty of about 10% on the measurements of $`v\mathrm{sin}i`$ . For a few stars of the sample, most of them presenting evidence of high rotation, the $`v\mathrm{sin}i`$ value was determined on the basis of the spectral synthesis carried out through this work.

### 2.2 Spectroscopic observations

The spectral region around the lithium line at 6707.81 Å was observed with two different telescopes. For northern stars, high-resolution spectra of the lithium region were acquired with the AURELIE spectrograph (Gillet at al., 1994) mounted at the 1.52 m telescope of the Observatoire de Haute Provence (France). The spectrograph used a cooled 2048-photodiode detector forming a 13 $`\mu `$m pixel linear array. A grating with 1800 lines/mm was used, giving a mean dispersion of 4.7 Å/mm and a resolving power around 45 000 (at 6707 Å). For this instrumentation and the selected set-up, the spectral coverage was about 120 Å. The signal to noise ratio was always better than 50. For southern stars, high-resolution spectra were acquired with the Coudé Echelle Spectrometer (CES) in the long camera mode (Kapper & Pasquini, 1996), mounted at the 1.44 m CAT telescope, at La Silla, ESO. A RCA high resolution CCD with 640 x 1024 pixels was used as detector, with a pixel size of 15 x 15 $`\mu `$m. The dispersion was around 1.9 Å/mm and the resolving power was about 95 000 (at 6707 Å). The spectral coverage for this instrument was about 70 Å and the signal to noise ratio was always better than 80. For both observing runs, thorium lamps were observed before and after each stellar observation for wavelength calibration, whereas 3 series of flat-field using an internal lamp of tungstene were obtained during each night of observation. Flat-field corrections and wavelength calibrations were performed using the MIDAS package.

## 3 Lithium abundances and spectrum synthesis

In the present study, we adopted the spectral analysis method used by Lèbre at al. (1999) and Jasniewicz et al. (1999), in order to derive the lithium abundances from the resonance LiI line (6707.81 Å). We refer to these authors for a description of our abundance analysis assuming LTE. Nevertheless a few important points are described below.

To compute the more appropriated synthetic spectra to fit the observations, a first guess for the stellar parameters of our list of giants has been estimated. Thus we have looked into the literature and found effective temperatures ($`T_{\mathrm{eff}}`$ ), gravities (log $`g`$ ) and often metallicites (\[Fe/H\]) for 37 stars (over 54), mainly in Cayrel de Strobel et al. (2001), in Allende Prieto and Lambert (1999), and in Alonso et al. (1999). For the rest of our sample (17 stars without any published stellar parameters determination), we first estimated $`T_{\mathrm{eff}}`$ from (B-V) colour index (Flower, 1996), and we have set log $`g`$ = 3.0, a microtuburlence velocity of 2 $`\mathrm{km}\mathrm{s}^1`$ and \[Fe/H\] = 0.0 which are values commonly adopted for Pop. I giant stars.

Synthetic spectra were then computed from these stellar parameters and from the same line list and model atmospheres used in Lèbre at al. (1999) and Jasniewicz et al. (1999). The stellar parameters of all the giants ( mainly $`T_{\mathrm{eff}}`$ , \[Fe/H\] and $`v\mathrm{sin}i`$ ) were then corrected, when necessary, in order to improve the fit quality of the several Fe I and other metallic lines found in the observed spectral range. The final adopted stellar parameters and the derived Li abundances are presented in Table 1 for all the stars of the sample, except HD 156015 which was too cool to accurately derive its effective temperature with the adopted procedure (however, the absence of Li signature is clear for this star). As already discussed by Lèbre et al. (1999), the major source of uncertainty for this abundance analysis is due to errors in the determination of the $`T_{\mathrm{eff}}`$ . We estimated that effective temperature were derived with an uncertainty smaller than $`\pm `$ 200 K. This leads to an error of less than 0.2 dex on the derived metallicities and lithium abundances. We also checked with the work of Carlsson et al. (1994) that non-LTE effects can be neglected for the studied stars since these effects are always much smaller than 0.1 dex.

## 4 Results and Discussion

### 4.1 Rotation and ALi in the Hertzsprung gap

The first step for the present analysis was the construction of the HR diagram to locate the evolutionary stage of the stars composing the working sample. In fact, such a procedure is important because in previous studies the only criterion for giant star classification was the spectral type. For this purpose, we have used HIPPARCOS (ESA 1997) trigonometric parallax measurements and $`V`$ magnitudes to compute absolute magnitudes and luminosities. Bolometric corrections $`BC`$ were determined from Flower (1996). Evolutionary tracks were computed from the Toulouse-Geneva code for stellar masses between 1 and 4 $`M_{}`$and solar metallicity (see also do Nascimento et al. 2000 for complementary informations). The HR diagrams with the referred evolutionary tracks are displayed in Figs. 1 and 2. In addition, these figures show the behaviors of the lithium abundance and of the rotational velocity $`v\mathrm{sin}i`$ , respectively. In these diagrams the dashed line indicates the evolutionary region where the subgiant branch starts, corresponding to the hydrogen exhaustion in stellar central regions, whereas the dotted line represents the beginning of the ascent along the red giant branch. Except for a few stars located along the subgiant branch, the large majority of stars in the present sample are effectively giants evolving prior to the ascent of the red giant branch.

Figure 1 shows the well established rotational discontinuity for giants of luminosity class III (e.g. De Medeiros and Mayor 1990). As already shown by these authors, giants blueward of the spectral G0III (corresponding to the location of the rotational discontinuity, around log($`T_{\mathrm{eff}}`$ ) $``$ 3.75), show a wide range of rotational velocity from a few $`\mathrm{km}\mathrm{s}^1`$ to about one hundred times the solar rotation. Giants redward of such spectral type (G0III) are essentially slow rotators, except for the synchronized binary systems and a dozen of late–G and K single giants for which the origin of their high rotation is not yet fully understood. At least for stars with turnoff mass greater than 1.5 $`M_{}`$, the sudden decline in rotation appears simply to reflect the rapid increase of the moment of inertia as the star evolves across this region of the HR diagram.

Figure 2 presents the behavior of lithium abundance, with a sudden decline in ALi near the same spectral region where the discontinuity in rotation is observed. Nevertheless, it is quite clear from Figs. 1 and 2 that the behavior of rotation and lithium content in the Hertzsprung gap, as well as the location of both discontinuities, strongly depends on stellar mass. Such a fact shows that, for giant stars, it is not correct to define the location of the discontinuities in rotation and in ALi at the same spectral type, as suggested by previous works based only on the analysis of the distributions of ALi and $`v\mathrm{sin}i`$ versus spectral type or (B$``$V) color index. Single stars located on the blue side of the Hertzsprung gap show a large spread in the lithium abundance, with values of ALi ranging from less than 0.5 to 3.0 dex. In contrast, stars located at the red side show essentially low values of ALi, reflecting the dilution effects along this spectral region.

Clearly, one observes from Fig. 2 that stars with 2 to 3 $`M_{}`$located redward of G0III (log($`T_{\mathrm{eff}}`$ ) $``$ 3.75) show a factor of dilution of at least 600, which is far in excess from the theoretical predictions. Such a fact points to an extra-mixing mechanism occuring before the beginning of the ascent of the red giant branch. The less massive giants show a similar disagreement between predicted and observed abundances. Brown et al. (1989) have measured lithium abundances for a large sample of late-G to K giants, with mean mass around 1.5 $`M_{}`$. For such stars, they found the mean lithium abundance lower than 0.0 dex, that is at least $``$ 1.5 dex below the value predicted by standard theory. In Fig. 2 of the present paper, by taking into consideration the abundance of lithium for stars with masses between 1.5 $`M_{}`$and 3.0 $`M_{}`$, located at the blue side of the Hertzsprung gap, and referring to a mean value of ALi lower than 0.0 dex as estimated by Brown et al. (1989), we observe that the factor of dilution for this interval of mass is even more important than that observed for stars with 1.5 $`M_{}`$.

It is interesting to point out that the binary systems seem to follow the same tendency presented by single stars. Redward of the spectral type G0III, the binary systems show essentially low values of ALi, reflecting also the effects of dilution along the Hertzsprung gap.

### 4.2 The connection lithium - rotation in the Hertzsprung gap

The dependence of lithium content on rotation along the Hertzsprung gap has been already reported in different works (Wallerstein et al., 1994; Alschuler 1975). For the present sample of stars, this trend is also clearly observed from Fig. 3, which displays lithium abundance as a function of rotational velocity. Stars with high rotation also present high values of ALi. In addition, the present data show a trend also observed for other luminosity classes, as well as for other spectral regions of giants. For slow rotators - in the present situation, for stars with $`v\mathrm{sin}i`$ lower than about 20 $`\mathrm{km}\mathrm{s}^1`$ \- there is a large spread of the values of ALi, with a dispersion of at least 3 magnitudes. Such a feature indicates that, also in the Hertzsprung gap, rotation is not the unique parameter driving Li abundance. Let us recall that, in spite of the limited number of binary systems with high $`v\mathrm{sin}i`$ value, the dependence of lithium on rotation seems to follow the same trend observed for the single stars. We have also analyzed the relation metallicity versus $`v\mathrm{sin}i`$ effects, looking for possible effects of metallicity on the spread observed in the ALiversus rotation relation, by using data presented in Table 1. No clear trend arises from such an analysis.

### 4.3 The connection lithium - deepening of the convective envelope

Because the level of dilution of lithium depends on the level of convection in the stars, we analyse here the behavior of lithium abundance as a function of the deepening of the convective zone for our working sample. Such an analysis sounds interesting because most of the stars are crossing the Hertzsprung gap, where the convective zone is predicted to reach its most important development. For our purposes, we have first estimated the mass of each star $`M_{}`$ from the HR diagram constructed for Figs. 1 and 2. Then, we have estimated the deepening in mass of the convective zone $`M_{CZ}`$, according to the recipe applied by do Nascimento et al. (2000) and do Nascimento et al. (2003). In short, with the mass of each star in hand, we have placed such mass in the convective zone mass deepening $`M_{\mathrm{cz}}/M_{}`$ versus Teff diagram from do Nascimento et al. (2000) and estimated the parameter $`M_{\mathrm{cz}}/M_{}`$ for each star. The individual values of $`M_{}`$ and $`M_{\mathrm{cz}}/M_{}`$ are listed in Table 1. The star HD 224342 is outside this analysis because of its very large and uncertain luminosity and mass. Figure 4 displays the deepening (in mass) of the convective envelope versus effective temperature, for single and binary stars. In this figure, the presence of the binary stars should be analysed with caution due to the large uncertainties in the estimation of mass for stars in binary sytems using evolutionary tracks. However it is clear from this figure that the sudden decline in ALi in the Hertzsprung gap is directly associated with the rapid increase of the convective envelope in such region. Most of the stars with high Li content present an undeveloped convective zone, whereas stars with low Li content show a developed convective zone.

## 5 Conclusions

The present study brings precise abundances of Li for giant stars of luminosity class III evolving across the Hertzsprung gap. By combining these data with rotational velocity and other stellar parameters such as luminosity and effective temperature, we show a trend for the existence of a sudden decline in ALi located around log($`T_{\mathrm{eff}}`$ ) $``$ 3.75, and corresponding to the spectral type G0III, where the well established rotational discontinuity for giants is defined. Blueward of this region, most of the stars with mass greater than about 1.5 $`M_{}`$present a trend for high Li content, whereas for stars redward of this region Li is essentially diluted. The scenario of Li content in single stars seems to be followed by binaries. The dependence of Li content on rotation, in the sense that fast rotators show high ALi, is also confirmed by the present set of data. In addition, the large spread in ALi for slow rotators observed for other luminosity classes and other spectral regions occurs also in the Hertzsprung gap. Metallicity seems to have no effect on such a spread. Finally, we have combined the ALi values with theoretical predictions on the convective zone deepening to analyse the extent of the level of convection on the Li dilution. It is clear from this analysis that stars with high Li content have mostly undeveloped convective zones, whereas stars with developed convective zone present low Li content.

###### Acknowledgements.

JRM warmly acknowledges colleagues of the Cote d’Azur Observatory at Nice, where a large part of this work was prepared during his stay in France. This work has been supported by continuous grants from the CNPq Brazilian Agency and by funds from the french Programme National de Physique Stellaire (INSU/CNRS). J.D.N.Jr. acknowledges the CNPq for the fellowship PROFIX 540461/01-6. We also gratefully acknowledge C. Ferrari for her wonderful software tiramisu.


format: text

text:
# 1 Two types of fiber arrangement for nonlinearity compensation with OPC: mirror-symmetric (left) and translationally symmetric (right).

Two Means of Compensating Fiber Nonlinearity

Using Optical Phase Conjugation

Haiqing Wei and David V. Plant

Department of Electrical and Computer Engineering
McGill University, Montreal, Canada H3A-2A6

hwei1@po-box.mcgill.ca

## Abstract

Two fiber lines may compensate each other for nonlinearity with the help of optical phase conjugation. The pair of fiber lines and the optical signals in them may be either mirror-symmetric or translationally symmetric about the conjugator.

© 2024 Optical Society of America
OCIS codes: (060.2330) Fiber optics communications; (190.4370) Nonlinear optics, fibers

The growing demand of higher capacity over longer transmission distances has seen fiber nonlinearity as a major limiting factor in modern optical transmission systems . Among the methods under investigation, nonlinearity compensation using optical phase conjugation (OPC) has emerged as an effective means of suppressing the nonlinear impairments . This paper shall discuss two types of fiber arrangement with respect to the OPC, as shown in Fig.1. In one type of arrangement, the fiber parameters and the signal intensity are in scaled mirror symmetry about the OPC. While the other type is characterized by scaled translational symmetry.

A mirror-symmetric link may consist of a fiber line on the left stretching from $`z=L/R`$ to $`z=0`$, $`L>0`$, $`R>0`$, followed by an OPC, then a fiber line on the right stretching from $`z=0`$ to $`z=L`$. The two fiber lines may carry wavelength-division multiplexed (WDM) signals $`_nA_n(z,t)\mathrm{exp}\left[i^z\beta (\zeta ,\omega _n)𝑑\zeta i\omega _nt\right]`$ and $`_nA_n^{}(z,t)\mathrm{exp}\left[i^zq(\zeta ,\omega _n^{})𝑑\zeta i\omega _n^{}t\right]`$ respectively, where $`n𝐙`$, $`\omega _n=\omega _0+n\mathrm{\Delta }`$ and $`\omega _n^{}=\omega _0^{}+n\mathrm{\Delta }`$ are the center frequencies of the WDM channels, $`\mathrm{\Delta }>0`$, $`\omega _0`$ is not necessarily equal to $`\omega _0^{}`$, $`A_n`$ and $`A_n^{}`$ are the slow-varying envelopes, while $`\beta (z,\omega )`$ and $`q(z,\omega )`$ are the $`z`$-dependent propagation constants. Being neglected is the random polarization-mode dispersion (PMD) effect. And for mathematical simplicity, all optical signals are assumed co-linearly polarized in the fibers. Define short-hand notations $`\beta ^{(k)}(z,\omega )\stackrel{\mathrm{def}}{=}\frac{^k}{\omega ^k}\beta (z,\omega )`$, and $`q^{(k)}(z,\omega )\stackrel{\mathrm{def}}{=}\frac{^k}{\omega ^k}q(z,\omega )`$, $`k>0`$. The dynamics of signal propagation in the two fiber lines is governed by two groups of coupled partial differential equations  respectively,

$`{\displaystyle \frac{A_n}{z}}i\beta _{2n}(z,i{\displaystyle \frac{}{t}})A_n+{\displaystyle \frac{\alpha _n(z)}{2}}A_n=i\gamma (z){\displaystyle \underset{k}{}}{\displaystyle \underset{l}{}}A_kA_lA_{k+ln}^{}\mathrm{exp}[i\theta _{kln}(z)]+{\displaystyle \underset{m}{}}g_{mn}(z)|A_m|^2A_n,`$ (1)
$`{\displaystyle \frac{A_n^{}}{z}}iq_{2n}(z,i{\displaystyle \frac{}{t}})A_n^{}+{\displaystyle \frac{\alpha _n^{}(z)}{2}}A_n^{}=i\gamma ^{}(z){\displaystyle \underset{k}{}}{\displaystyle \underset{l}{}}A_k^{}A_l^{}A_{k+ln}^{}\mathrm{exp}[i\theta _{kln}^{}(z)]+{\displaystyle \underset{m}{}}g_{mn}^{}(z)|A_m^{}|^2A_n^{},`$ (2)
$`n𝐙`$, where for the first fiber line, $`\alpha _n`$ is the attenuation coefficient around $`\omega _n`$, $`\gamma `$ is the Kerr nonlinear coefficient, $`g_{mn}`$ is the Raman coupling coefficient from the $`m`$th to the $`n`$th channels, $`\theta _{kln}(z)\stackrel{\mathrm{def}}{=}^z[\beta (\zeta ,\omega _k)+\beta (\zeta ,\omega _l)\beta (\zeta ,\omega _{k+ln})\beta (\zeta ,\omega _n)]𝑑\zeta `$ is the phase mismatch among the mixing waves, and the functional operator $`\beta _{2n}`$ is defined as $`\beta _{2n}(z,i\frac{}{t})\stackrel{\mathrm{def}}{=}_{k=1}^+\mathrm{}\frac{1}{k!}\beta ^{(k)}(z,\omega _n)\left(i\frac{}{t}\right)^k\beta ^{(1)}(z,\omega _0)\left(i\frac{}{t}\right)`$. The parameters $`\alpha ^{}`$, $`\gamma ^{}`$, $`g_{mn}^{}`$, $`\theta _{kln}^{}`$ and the operator $`q_{2n}(z,i\frac{}{t})`$ are similarly defined for the second fiber line. It is an easy exercise to show that equations (1) reduce to (2), when the parameters satisfy the following rules of correspondence,

$`\beta ^{(2)}(z,\omega _0+\omega )`$ $`=`$ $`Rq^{(2)}(Rz,\omega _0^{}\omega ),\omega 𝐑,`$ (3)
$`\alpha _n(z)`$ $`=`$ $`R\alpha _n^{}(Rz),n𝐙,`$ (4)
$`\gamma (z)`$ $`=`$ $`R\gamma ^{}(Rz)|C|^2,`$ (5)
$`g_{mn}(z)`$ $`=`$ $`Rg_{m,n}^{}(Rz)|C|^2,m,n𝐙,`$ (6)
$`z[0,L/R]`$, $`C0`$ being a constant, and the envelope functions are related as $`A_n(z,t)=CA_n^{}(Rz,t)`$, $`n𝐙`$. Physically, it says that the two fiber lines compensate each other for dispersion and nonlinearity. Optical signals $`A_n(L/R,t)`$, $`n𝐙`$, entering the first fiber line may be dispersed and nonlinearly distorted to become $`A_n(0,t)`$, $`n𝐙`$, which are converted into $`A_n^{}(0,t)=A_n^{}(0,t)/C^{}`$, $`n𝐙`$, by the OPC. The second fiber line will then propagate the optical signals in a reversed manner with respect to the first. The final outputs signals $`A_n^{}(L,t)=A_n^{}(L/R,t)/C^{}`$, $`n𝐙`$, are exact replicas of the initial signals up to complex conjugation. It is noted that parts of one fiber line would amplify light in correspondence to the attenuation in parts of the other, and vice versa. A specialty fiber may be chosen with parameters satisfying equations (3,5,6) to be the scaled mirror image of a transmission fiber which usually attenuates light. At the same time, erbium doping or Raman pumping should be employed to obtain the gain specified by (4).

A link with translational symmetry could be constructed to cancel weak nonlinearities up to the first order perturbation. Consider two fiber lines with opposite Kerr and Raman nonlinear coefficients but identical linear parameters. If (1) with $`z[L,0]`$ describe the signal propagation in one fiber line, then the signal dynamics in the other would be governed by similar equations with negative $`\gamma `$ and $`g`$ coefficients,

$`{\displaystyle \frac{B_n}{z}}i\beta _{2n}(zL,i{\displaystyle \frac{}{t}})B_n+{\displaystyle \frac{\alpha _n(zL)}{2}}B_n=i\gamma (zL){\displaystyle \underset{k}{}}{\displaystyle \underset{l}{}}B_kB_lB_{k+ln}^{}\mathrm{exp}[i\theta _{kln}(zL)]`$
$`{\displaystyle \underset{m}{}}g_{mn}(zL)|B_m|^2B_n,0zL,`$ (7)
which take the input $`B_n(0,t)`$, $`n𝐙`$ and give the output $`B_n(L,t)`$, $`n𝐙`$. When the signal intensity is not very high, so that the nonlinearity is weak and treated with perturbation theory, the output from each fiber line is a linearly dispersed version of the input, plus nonlinear distortions expanded in power series of the $`\gamma `$ and $`g`$ coefficients. By neglecting the higher order powers and keeping only the terms linear in $`\gamma `$ or $`g`$, it can be seen that the two fiber lines induce opposite nonlinear distortions to otherwise the same, linearly dispersed signals. If the overall dispersion of each line is compensated to zero and the signal loss is made up by a linear optical amplifier, then the two lines in cascade would comprise a transmission line with fiber nonlinearity annihilated up to the first order perturbation. The problem is that an optical fiber with negative nonlinear coefficients does not exist naturally. Fortunately, it can be simulated by a regular fiber with the help of OPC. Take a regular fiber with parameters $`(q,\alpha ^{},\gamma ^{},g^{})`$ that satisfy,

$`q^{(2)}(z,\omega _0^{}+\omega )`$ $`=`$ $`R\beta ^{(2)}(RzL,\omega _0\omega ),\omega 𝐑,`$ (8)
$`\alpha _n^{}(z)`$ $`=`$ $`R\alpha _n(RzL),n𝐙,`$ (9)
$`\gamma ^{}(z)`$ $`=`$ $`R\gamma (RzL)|C|^2,`$ (10)
$`g_{mn}^{}(z)`$ $`=`$ $`Rg_{m,n}(RzL)|C|^2,m,n𝐙,`$ (11)
$`z[0,L/R]`$. The signal propagation in the regular fiber is then governed by,

$$\frac{B_n^{}}{z}iq_{2n}(z,i\frac{}{t})B_n^{}+\frac{\alpha _n^{}(z)}{2}B_n^{}=i\gamma ^{}(z)\underset{k}{}\underset{l}{}B_k^{}B_l^{}B_{k+ln}^{}\mathrm{exp}[i\theta _{kln}^{}(z)]+\underset{m}{}g_{mn}^{}(z)|B_m^{}|^2B_n^{},$$
(12)
$`z[0,L/R]`$, which are solved by $`B_n^{}(z,t)=CB_n^{}(Rz,t)`$, $`n𝐙`$, and turn the input $`B_n^{}(0,t)=CB_n^{}(0,t)`$, $`n𝐙`$, into the output $`B_n^{}(L/R,t)=CB_n^{}(L,t)`$, $`n𝐙`$. The regular fiber equipped with OPC at its two ends takes the input $`B_n(0,t)`$, $`n𝐙`$ and gives the output $`B_n(0,t)`$, $`n𝐙`$. That fulfils the function of the fictitious fiber with negative nonlinearity. The OPC at the output end of the regular fiber may be omitted in practice, as most applications would not differentiate between a signal and its conjugate.

It is noted that each fiber line on one side of the OPC is not necessarily one fiber span, and the signal intensity does not have to evolve monotonically either. Both methods work fine when each side of the OPC consists of multiple fiber spans with optical amplifiers boosting the signal power, although the added noise makes perfect nonlinearity compensation impossible. Using a commercial software, computer simulations have been carried out to test the proposed methods of nonlinearity compensation. For the mirror setup, the test link consists of a specialty fiber, an OPC, and a transmission fiber $`200`$ km long, with loss $`\alpha =0.2`$ dB/km, dispersion $`D=8`$ ps/nm/km, dispersion slope $`S=0.08`$ ps/nm2/km, effective mode area $`A_{\mathrm{eff}}=50`$ $`\mu `$m2, Kerr and Raman coefficients that are typical of silica glass. The specialty fiber is made of the same material, but with parameters $`(\alpha ^{},D^{},S^{})=20\times (\alpha ,D,S)`$ and $`A_{\mathrm{eff}}^{}=12.5`$ $`\mu `$m2. The nonlinearity of the specialty fiber can be switched on and off. Amplifier noise is added at the two ends of the link. The input are four WDM channels spaced by $`100`$ GHz, return-to-zero modulated at $`10`$ Gb/s with $`33\%`$ duty. The pulses peak at $`100`$ mW when entering the transmission fiber. For a WDM system with the span loss so large and the input optical power so high, the output signals would be distorted heavily and become unusable, were there no nonlinearity compensation . By contrast, the transmission system becomes virtually penalty-free with our scheme of mirror-symmetric nonlinearity compensation . The test system in translational symmetry is constructed with ten $`100`$-km spans on one side of the OPC using the same transmission fiber as in the mirror setup. Each span is ended by an erbium-doped fiber amplifier (EDFA) with $`20`$ dB gain, noise figure $`4`$ dB, and a dispersion compensating module (DCM) with negligible nonlinearity. The DCM perfectly compensates the $`D`$ and $`S`$ of the fiber span. On the other side of the OPC are ten spans of transmission fiber with parameters $`(\alpha ^{},D^{},S^{},A_{\mathrm{eff}}^{},\gamma ^{},g^{})=(\alpha ,D,S,A_{\mathrm{eff}},\gamma ,g)`$. The loss and dispersion of each span are also fully compensated by an EDFA and a DCM. The EDFA noise figure is still $`4`$ dB. The input RZ pulses peak at 10 mW when entering the transmission fiber. Firstly, the OPC is absent and the transmission result of the 20-span link is shown in the left-side graph of Fig.2. When the OPC is put back, the other graph in Fig.2 clearly demonstrates the effect of nonlinearity compensation.


text:
# Near-Infrared, Adaptive Optics Observations of the T Tauri Multiple-Star System

## 1 Introduction

Stars frequently form in small, gravitationally-bound groups of three or more (Ghez, Neugebauer, & Matthews, 1993; Mathieu, 1994), usually situated so that each perturbs the motion of the others significantly. In such systems the stellar orbits develop chaotically. Over time, gravitational interactions among the members and with the residue of the disks from which they formed, result in large orbital changes and even ejection, especially for the lighter components (see, e.g., Marchal, 1990). As main-sequence stars are always seen in configurations that are quite different from those of young multiple stars, direct observation of orbital evolution in young multiple-star systems has long been sought.

Thus the report by Loinard, Rodríguez, & Rodríguez (2003) was greeted with wide interest: one of the components of T Tauri has, in the last few years, undergone a dramatic change in orbital velocity, possibly resulting in ejection from the system. T Tau has long been the archetype of young solar-analog stars, and more recently has become an archetype of young multiple star systems. Since the discovery two decades ago that T Tau is not single, but has a companion 0 $`\stackrel{}{\mathrm{.}}`$7 to the south (Dyck, Simon, & Zuckerman, 1982), high-resolution infrared observations using speckle interferometry or adaptive optics have been used to resolve the system into three stars and to measure accurately the motions arising from their orbits: the visible, classic T Tau N, and two heavily-extinguished companions, T Tau Sa and T Tau Sb, separated by about 0 $`\stackrel{}{\mathrm{.}}`$1 (Koresko, 2000; Köhler, Kasper, & Herbst, 2000; Duchêne, Ghez, & McCabe, 2002). Meanwhile, radio astronomers used the VLA at 2 cm wavelength to resolve the system also into two compact components, one coincident with T Tau N and a brighter one close to T Tau S (Schwartz, Simon, & Campbell, 1986). The southern radio component exhibited motion that until recently was consistent with orbital motion around T Tau Sa (which does not appear in the radio images) and has been identified with T Tau Sb (Johnston et al., 2003; Loinard, Rodríguez, & Rodríguez, 2003; Smith et al., 2003). The radio source is unresolved in VLBI images (diameter $`<`$ 0.5 mas; Smith et al., 2003), variable in brightness, and exhibits strong and variable polarization in its radio emission (Smith et al., 2003; Johnston et al., 2003). This is the object that has apparently been ejected.

Here we present new near-infrared observations which we use to show that the southern radio component cannot be the same as the star T Tau Sb. We suggest that T Tau has been a quadruple stellar system, with the lowest-mass member probably ejected during the past few years.

## 2 Observations

We observed T Tau under good observing conditions on 2002 December 24 with the Palomar Observatory Hale 200-inch telescope, the Palomar Adaptive Optics (PALAO) high-order image-correction system (Troy et al., 2000), and the Palomar High-Angular-Resolution Observer (PHARO) near-infrared camera/spectrometer (Hayward et al., 2001). The brightest visible component of the system, T Tau N (V = 9.6), served as the phase reference for the AO system. We also observed SAO 76481 as flux calibrator and point-spread-function (PSF) standard immediately after taking our images of the T Tau system. We chose the image scale to be 25 mas per pixel.

In order to calibrate the plate scale and orientation of the array, we used our observations of the multiple systems RW Aur (two components) and UX Tau (four components) and compared the measured separations and position angles with the corresponding measurements in White & Ghez (2001). We chose multiple systems whose components have separations wide enough that orbital motion should be negligible in the time period between the observations by White & Ghez (2001) (1996-1997) and ours (2002.98). The plate scale of our images comes out to 24.7 $`\pm `$ 0.4 mas per pixel, and the orientation of the detector array is such that north is up and east is to the left with $`\pm `$ 1$`\mathrm{°}`$ accuracy. Thus position angles inferred from our data have uncertainty $`\pm `$ 1$`\mathrm{°}`$.

Figure 1 is the image of T Tau which we took in a narrow-band filter ($`\mathrm{\Delta }\lambda `$ = 0.06 $`\mu `$m) in the infrared K band centered at 2.26 $`\mu `$m. (The components of T Tau are bright enough at near-infrared wavelengths that the use of the usual broadband filters is inconvenient, as will be seen below). A total of 20 frames, each with exposure time 1.8 seconds, were added to create this image. The resulting PSF is close to diffraction-limited, 102 mas in diameter (FWHM). The limiting 2.26 $`\mu `$m magnitude for parts of the image well separated from the brighter stars is 14.7 (5 $`\sigma `$ per pixel). Within 200 mas (two PSF widths) of the stars T Tau Sa and Sb, diffraction rings reduce the limiting magnitude to a range from 10.2 to 11.2. All three of the previously-identified stellar components of T Tau appear in Figure 1, well resolved from one another. Table 1 is a list of their positions and magnitudes. T Tau Sb lies 0 $`\stackrel{}{\mathrm{.}}`$107, at position angle 289$`\mathrm{°}`$, from T Tau Sa, and the magnitudes of the two southern components are about the same, with T Tau Sb brighter by a factor of 1.05. Since in late 2000 the corresponding separation and position angles amounted to 0 $`\stackrel{}{\mathrm{.}}`$092 and 267$`\mathrm{°}`$, respectively (Duchêne, Ghez, & McCabe, 2002), we probably see orbital motion of T Tau Sb around Sa (see also Figure 3).

The positions and fluxes of the three observed infrared components of the T Tau system were determined by applying a PSF fit to the data. Standard IDL procedures were used to perform a Gaussian fit to the given PSF, the one of SAO 76481, and to apply it to the target stars. The position uncertainty resulting from this method is estimated to be 0.2 pixels, which corresponds to 5 mas. The flux calibration was carried out by using the 2MASS K-band flux density of SAO 76481; we estimate the uncertainty of the absolute photometry to be $`\pm `$ 20 %, and that of the relative photometry to be a few percent.

Figure 2 is our Ks image ($`\lambda `$ = 2.145 $`\mu `$m, $`\mathrm{\Delta }\lambda `$ = 0.31 $`\mu `$m) of T Tau Sa and Sb, taken 5 minutes before our narrow-band K data. It is also a sum of 20 frames of 1.8 seconds exposure time each. Since T Tau N and our calibrator, SAO 76481, saturated even with this short integration time, we used the image of SAO 76481 observed with the 2.26 $`\mu `$m filter to perform PSF-fitting on this image. The position difference between Sa and Sb is nearly the same in this image and the one taken in the 2.26 $`\mu `$m filter; the offsets listed in Table 1 are averages of the two. No absolute flux calibration was possible, but relative photometry was: T Tau Sb is brighter than T Tau Sa by the factor 1.36, somewhat higher than we obtained with the narrow-band data. In another image taken about 30 minutes later with the Ks filter and an occulting mask over T Tau N, the flux ratio of T Tau Sa and Sb was the same as in the previous Ks image. The difference between this result and the flux ratio seen in the narrow-band K filter (Sb/Sa = 1.05) could be the result of differences in spectral features and degree of extinction between the two components.

Comparing our narrow-band K data to the results from late 2000 by Duchêne, Ghez, & McCabe (2002), T Tau Sa is fainter and T Tau Sb is brighter, each by about a magnitude at K. That T Tau Sb can be as bright as T Tau Sa in the K band means that the orbital motion of T Tau Sa relative to T Tau N is confused in all pre-1997 measurements because of the presence of T Tau Sb.

## 3 Analysis and Conclusions

Figure 3 is a plot of the radio/infrared positions of the southern components of T Tau between 1983 and 2003, in a frame of reference in which T Tau Sa is at rest. The radio positions in this plot differ very slightly from those determined by Loinard, Rodríguez, & Rodríguez (2003). They are based on the radio/infrared registration of T Tau N and include the (very small) orbital motion of T Tau Sa with respect to T Tau N; the uncertainty in the distance between origins of the radio and infrared reference frames is about 15 mas (radial, RMS). The relative positions of T Tau Sa and Sb were measured directly in the infrared images. Our new position for T Tau Sb is quite inconsistent with anything along the track of the radio source between 1998 and 2001: it misses by 50-100 mas, much larger than the combined uncertainty of about 16 mas. For example, the linear extrapolation of the two most recent radio-source positions, the point indicated by a empty diamond in Figure 3, lies 78 mas from the 2002.98 observed position of the infrared star T Tau Sb. If a single object were to have followed the path described by both the radio and infrared positions, its velocity in the plane of the sky would have had to undergo major changes twice, once between 1995 and 1998, and once after 2001. Moreover, in at least one of the cases, the single object would have had no obvious partner with which to interact. Thus the probability is much greater that the radio and infrared observations of the southern components of T Tau detect different objects. In anticipation of its detection in future infrared observations, we will refer to the radio source as T Tau Sc henceforth.

Also appearing in Figure 3 is a model orbit, resulting from a minimum-$`\chi ^2`$ fit to the projected trajectory and transverse velocities of T Tau Sc measured at the VLA between 1983 and 1996 (Loinard, Rodríguez, & Rodríguez, 2003; Johnston et al., 2003). In this exercise, we assumed T Tau Sa to be at one focus of the orbit. The fit to both trajectory and velocity improved after shifting all the VLA positions of T Tau Sc with respect to T Tau Sa by 4.5 mas to the east, an amount small compared to the uncertainty of the radio-infrared coordinate-system correlation. This small shift is present in the positions plotted in Figure 3. The agreement between model and observations is excellent; we suggest that the 4.5 mas shift is an improvement upon the previous determination of the infrared-radio reference-frame alignment. The parameters of the model orbit are as follows: period 19.1 years, semimajor axis 9.1 AU, assuming a distance of 140 pc (leading to a total mass of $`2.1M_{}`$ for T Tau Sa and T Tau Sc), axis inclination 59 degrees from the line of sight, eccentricity 0.57, argument of periastron 168.5 degrees, position angle 0.5 degrees for the line of nodes, and periastron epoch 1997.39. T Tau Sb is a young M1 $`\pm `$ 1 star (Duchêne, Ghez, & McCabe, 2002), for which the mass is much smaller than $`2.1M_{}`$; thus we have neglected it in this model.

What is the nature of T Tau Sc? Its small size and large (kilogauss) magnetic field, indicated by the VLBI observations (Smith et al., 2003), are consistent with a star-like object. We did not detect any other object besides T Tau N, Sa, and Sb in our near-infrared images. At the position extrapolated for T Tau Sc at epoch 2002.98 (see Figure 3), our upper limit for the detection of a point source is about 10.2 magnitudes at 2.26 $`\mu `$m. Thus it is either lighter than the eighth-magnitude T Tau Sb, or more heavily-extinguished, or both. The velocities of T Tau Sb and Sc within the T Tau system are larger than those of T Tau Sa, which, if they have been gravitationally bound, indicates that both objects are less massive than T Tau Sa. Therefore T Tau Sc is probably an ordinary young star or brown dwarf, similar in mass or lighter than T Tau Sb, though quite different from this star in its magnetic field strength.

The presence of two distinct objects, T Tau Sb and the radio source T Tau Sc, could provide a natural explanation for changes in each other’s state of motion. As Loinard, Rodríguez, & Rodríguez (2003) point out, the path of T Tau Sc obeys Kepler’s laws for an orbit around T Tau Sa from discovery through 1995, but departs strongly thereafter. The path of the infrared star T Tau Sb with respect to Sa shows a wide variation in area per unit time; between February 2000 and November 2000 the rate of area swept out is more than a factor of two greater than from late 1997 to early 2000, or from November 2000 to December 2002. The curvature and orientation of this path probably indicate that T Tau Sb has been, and is still, bound to T Tau Sa, but longer-term observations will be necessary to demonstrate this. In any case, both T Tau Sb and Sc experienced large orbital-motion changes in the 1998-2001 period, when they were closest together in projection. We suggest that their mutual gravitational interaction at this time has resulted in these orbit changes111If 20 years is a typical orbital period over the history of the T Tau S system, then there have been of the order of 105 revolutions since formation, during which there could have been only a few encounters that resulted in large orbital changes. It would seem that astronomers have been extremely fortunate to witness this event.. This would amend the suggestion by Loinard, Rodríguez, & Rodríguez (2003) that T Tau Sa has a close companion (separation $`<`$ 2 AU) whose interaction with the radio source T Tau Sc (which they identify as T Tau Sb) has resulted in the ejection of the latter; instead, T Tau Sb is a somewhat more distant companion that has probably caused the ejection of T Tau Sc, which is a fourth stellar or sub-stellar member in the T Tau system.

We are grateful to our telescope operator, Karl Dunscombe, and the instrument engineers, Rick Burruss and Jeff Hickey, for their outstanding support and help. This work is supported in part by NASA, through grants for the development and operation of the Space Infrared Telescope Facility Infrared Spectrograph: Jet Propulsion Laboratory (JPL) contract 960803 to Cornell University, and Cornell subcontracts 31419-5714 and 31419-5715 to the University of Rochester. L. D. K. and T. L. H. acknowledge support from the University Space Research Foundation under grant USRA 8500-98-014. This research has made use of the NASA/ IPAC Infrared Science Archive operated by JPL, California Institute of Technology (Caltech), under contract with NASA. It has also made use of the SIMBAD and VizieR databases, operated at CDS (Strasbourg, France), NASA’s Astrophysics Data System Abstract Service, and of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and IPAC/Caltech, funded by NASA and the National Science Foundation.


text:
# Geometry of large-scale low-energy excitations in the one-dimensional Ising spin glass with power-law interactions

## I Introduction

There have been several numerical attempts at finite temperatureKrzakala and Martin (2000); Palassini and Young (2000); Marinari and Parisi (2000); Katzgraber et al. (2001); Katzgraber and Young (2001, 2002); Houdayer et al. (2000); Houdayer and Martin (2000) and zero temperatureBray and Moore (1984); McMillan (1984a, b); Rieger et al. (1996); Hartmann (1999); Palassini and Young (1999); Hartmann and Young (2001); Carter et al. (2002) to better understand the nature of the spin-glass state for short-range spin glasses. These results are generally interpreted in terms of the two main theories for the spin-glass phase: the replica symmetry-breaking (RSB) picture,Parisi (1979, 1980, 1983); Mézard et al. (1987) and the “droplet picture”Fisher and Huse (1986, 1987, 1988); Bray and Moore (1986) (DP). RSB predicts that excitations involving a finite fraction of the spins cost only a finite energy in the thermodynamic limit, and that the fractal dimension of the surface of these excitations $`d_s`$ is equal to the space dimension $`d`$. This is in contrast to DP where a low-energy excitation (droplet) has an energy proportional to $`\mathrm{}^\theta `$, where $`\mathrm{}`$ is the characteristic length scale of the droplet and $`\theta `$ is a positive stiffness exponent. In addition, the surface of these excitations is fractal with $`d_s<d`$. More recently Krzakala and Martin,Krzakala and Martin (2000) as well as Palassini and Young,Palassini and Young (2000) suggest an intermediate picture (called “TNT” for trivial-nontrivial) in which droplets have a fractal surface with $`d_s<d`$, and their energy is finite in the thermodynamic limit. Which of the above pictures describes the spin-glass state correctly is still widely debated.

The RSB and TNT pictures require two stiffness exponents for the energy of large-scale excitations. There is convincing numerical evidence that changing the boundary conditions (e.g., from periodic to antiperiodic), which induces a domain wall, costs an energy which increases as $`\mathrm{}^\theta `$ with $`\theta >0`$. On the other hand, in the RSB and TNT pictures, the energy of droplets, created by thermal noise or by applying a perturbation for a fixed set of boundary conditions, varies as $`\mathrm{}^\theta ^{}`$ with $`\theta ^{}=0`$. By contrast, the DP makes the reasonable ansatz that $`\theta ^{}=\theta (>0)`$.

In a previous publication,Katzgraber and Young (2003) we studied the one-dimensional long-range Ising spin glass with power-law interactions. The model’s advantage is that large system sizes can be studied, in contrast to the short-range spin-glass models commonly used. The results of Ref. Katzgraber and Young, 2003 showed that the stiffness exponent $`\theta `$ for zero-temperature domain-wall excitations is positive and in fair agreement with analytical predictions from the droplet model.Fisher and Huse (1988); Bray et al. (1986); Katzgraber and Young (2003) However, the stiffness exponent $`\theta ^{}`$ for thermally induced droplet excitations is different and consistent with zero. Hence, at least for the range of system sizes studied, $`L512`$, the data of Ref. Katzgraber and Young, 2003 are consistent with both the TNT and RSB scenarios since they have $`\theta >\theta ^{}=0`$.

The purpose of the present paper is to estimate $`d_s`$, because this distinguishes between the RSB and TNT scenarios, since $`d_s=d`$ in RSB while $`d_s<d`$ in TNT. For short-range models, a droplet excitation forms a single connected piece, and so $`d_s`$ has to be zero in $`d=1`$. However, for long-range interactions, a droplet may consist of disconnected pieces,Fisher and Huse (1988) so a nontrivial value of $`d_s`$ is possible in $`d=1`$. We perform both finite-temperature Monte Carlo simulations and ground-state studies. Our results suggest that droplets are possibly compact in agreement with RSB, although the data are also consistent with a very small value of $`dd_s`$, which would be consistent with TNT.

In Sec. II we introduce the model, observables, and details of the Monte Carlo technique. Results at zero temperature are presented in Sec. III, and those at finite temperature are presented in Sec. IV. Our conclusions are summarized in Sec. V.

## II Model and Numerical Method

The Hamiltonian of the one-dimensional long-range Ising spin chain with power-law interactions is given by

$$=\underset{i,j}{}J_{ij}S_iS_j,$$
(1)
where the Ising spins $`S_i=\pm 1`$ are evenly distributed on a circular ring of length $`L`$ to ensure periodic boundary conditions. The sum is over all pairs of spins on the chain and the couplings $`J_{ij}`$ are given by

$$J_{ij}=c(\sigma )\frac{ϵ_{ij}}{r_{ij}^\sigma },$$
(2)
whereKatzgraber and Young (2003)

$$r_{ij}=\frac{L}{\pi }\mathrm{sin}\left(\frac{\pi |ij|}{L}\right)$$
(3)
is the straight-line distance between sites $`i`$ and $`j`$. The random part of the interactions $`ϵ_{ij}`$ is chosen according to a Gaussian distribution with zero mean and standard deviation unity, and the constant $`c(\sigma )`$ in Eq. (2) is chosenKatzgraber and Young (2003) to give a mean-field transition temperature $`T_c^{\mathrm{MF}}=1`$.

The one-dimensional long-range Ising spin chain has a very rich phase diagramFisher and Huse (1988); Bray et al. (1986); Katzgraber and Young (2003) in the $`d`$-$`\sigma `$ plane. Spin-glass behavior is controlled by the long-range part of the interaction if $`\sigma `$ is sufficiently small, and by the short-range part if $`\sigma `$ is sufficiently large. In this work we focus on the long-range behavior at $`\sigma =0.75`$ for whichKatzgraber and Young (2003) $`T_c>0`$ and the critical exponents are non-mean-field like. Using the exact relationBray et al. (1986); Fisher and Huse (1988) $`\theta =d\sigma `$, we expect $`\theta =0.25`$ for $`d=1`$, which is in moderate agreement with numerical resultsKatzgraber and Young (2003) for domain walls induced by a change in boundary conditions at $`T=0`$. By studying thermally induced droplet excitations Ref. Katzgraber and Young, 2003 also estimated $`\theta ^{}0`$, consistent with RSB and TNT.

In order to excite droplets at zero temperature we use the coupling-dependent ground-state perturbation method described elsewhere.Palassini and Young (2000); Palassini et al. (2002); Hartmann and Young (2002) First, we compute the ground-state configuration $`\{S_i^{(0)}\}`$. Then we perturb the couplings $`J_{ij}`$ by the following amount:

$$\mathrm{\Delta }(ϵ)=\frac{2ϵ}{N}\underset{i,j}{}\frac{[J_{ij}^2]_{\mathrm{av}}}{(T_c^{MF})^2}S_i^{(0)}S_j^{(0)}S_iS_j,$$
(4)
where $`ϵ`$ is a coupling constant and $`[\mathrm{}]_{\mathrm{av}}`$ represents a disorder average. The (total) energy of the unperturbed ground state then increases by exactly $`ϵ`$, whereas the energy of any other state $`\alpha `$ will increase by the lesser amount of $`ϵq_l`$, where $`q_l`$ is the link overlap between the unperturbed ground state and a state $`\alpha `$:

$$q_l=\frac{2}{N}\underset{i,j}{}\frac{[J_{ij}^2]_{\mathrm{av}}}{(T_c^{MF})^2}S_i^{(\alpha )}S_j^{(\alpha )}S_i^{(0)}S_j^{(0)}.$$
(5)
In previous work $`q_l`$ has been defined for nearest-neighbor models in which the sum is over nearest-neighbor pairs. Here we have generalized the link overlap to long-range models in a natural way. Because the coupling constant $`ϵ`$ is of order unity and not of order $`L`$ only low-energy excitations can be generated. We compute the new ground state of the perturbed system and record the link overlap between the old and new ground states. In the context of zero-temperature simulations the term “link overlap” will hereafter refer to the link overlap between the perturbed and unperturbed ground states.

Ground states are calculated using the parallel tempering Monte Carlo methodHukushima and Nemoto (1996); Marinari (1998) (at very low temperatures) as described in Refs. Katzgraber and Young, 2003 and Moreno et al., 2003. The parameters used in the $`T=0`$ simulations are shown in Table 1. For each value of $`L`$ and $`ϵ`$ we compute $`10^4`$ disorder realizations. We find that for $`\sigma =0.75`$, when the model is in the long-range phase, the efficiency of the used algorithm to calculate ground states scales as $`L^z`$ with $`z=2.9\pm 0.3`$. For the current project this translates to a total CPU time of 70 CPU years. Curiously, for $`\sigma =2.50`$, for which the interactions are effectively short range so frustration is minimal in the $`d=1`$ model studied here, the algorithm performs poorly with the equilibration time varying as $`\mathrm{exp}(aL)`$, with $`a=0.13\pm 0.02`$. It would be useful to understand intuitively the reasons for this.

One quantity that we study at $`T=0`$ is the link overlap, averaged over all samples. To see how this varies with sizePalassini and Young (2000) consider a large cluster of excited spins. This has a characteristic energy of order $`L^\theta ^{}`$, which is to be compared with the energy gained from the perturbation $`ϵ(1q_l)ϵL^{(dd_s)}`$. There is a distribution of cluster energies which we assume to have a finite weight at the origin, so the probability that the perturbation will create the excitation is $`ϵL^{(\theta ^{}+dd_s)}`$. When this occurs $`1q_lL^{(dd_s)}`$, and so on averagePalassini and Young (2000); Palassini et al. (2002); Hartmann and Young (2002)

$$[1q_l]_{\mathrm{av}}=ϵL^{\mu _l}(a+bL^c),$$
(6)
where

$$\mu _l=\theta ^{}+2(dd_s)$$
(7)
and we have added a correction to scaling term $`bL^c`$.

In RSB we have $`\mu _l=0`$ so Eq. (6) tends to a constant for $`L\mathrm{}`$, whereas in DP and TNT $`\mu _l>0`$ so Eq. (6) tends to zero in this limit.

In addition, we consider averages over only those samples in which a large excitation is generated,Hartmann and Young (2002); Palassini et al. (2002) comprising a finite fraction of spins. The criterion we take is $`|q|0.50`$. Averaging just over these samples givesHartmann and Young (2002); Palassini et al. (2002)

$$[1q_l]_{\mathrm{av}}^{}=L^{(dd_s)}(a+bL^c),$$
(8)
the prime representing the restricted average. Equation (8) follows trivially from the arguments presented in the derivation of Eq. (6) with the probability factor $`ϵL^{(\theta ^{}+dd_s)}`$ replaced by unity. We expect that $`[1q_l]_{\mathrm{av}}^{}`$ will be independent of $`ϵ`$.

In order to study droplet geometries at finite temperatures, we compute the distribution of the link overlap $`q_l`$ between two replicas $`\alpha `$ and $`\beta `$ of the system with the same disorder:

$$q_l=\frac{2}{N}\underset{i,j}{}\frac{[J_{ij}^2]_{\mathrm{av}}}{(T_c^{MF})^2}[S_i^{(\alpha )}S_j^{(\alpha )}S_i^{(\beta )}S_j^{(\beta )}]_{\mathrm{av}}.$$
(9)
Here $`\mathrm{}`$ represents a thermal average, and $`[\mathrm{}]_{\mathrm{av}}`$ represents a disorder average. From the finite-size scaling argumentsKatzgraber et al. (2001) used to derive Eq. (6) we expect that the variance of the distribution of the link overlap scales as

$$\mathrm{Var}(q_l)=L^{\mu _l}\left(a+bL^c\right).$$
(10)
Note that in RSB $`\mu _l=0`$ so $`\mathrm{Var}(q_l)`$ tends to a constant for $`L\mathrm{}`$. However, $`\mu _l>0`$ in DP (since $`\theta ^{}=\theta >0`$ and $`d_s<d`$) and in TNT (since $`d_s<d`$). The $`bL^c`$ term is a correction to scaling, which turns out to be necessary since the data cannot be fitted without it.

To speed up equilibration of the finite-$`T`$ simulations we use the parallel tempering Monte Carlo method.Hukushima and Nemoto (1996); Marinari (1998) We test for equilibration using the criterion developed earlier,Katzgraber et al. (2001) now generalizedKatzgraber and Young (2003) for the Hamiltonian in Eq. (1). For all sizes, the lowest temperature used is $`T=0.05`$, well below $`T_c0.63`$.Katzgraber and Young (2003); Leuzzi (1999) The highest temperature is $`1.70`$ which is well above the mean-field critical temperature ( $`T_c^{\mathrm{MF}}=1`$) and so the spins equilibrate fast there. We choose the spacing between the temperatures such that the acceptance ratios for the global moves are around $`0.30`$. Parameters of the finite-$`T`$ simulations are summarized in Table 2.

To summarize, for $`L\mathrm{}`$ all the quantities that we calculate, $`[1q_l]_{\mathrm{av}}`$ in Eq. (6), $`[1q_l]_{\mathrm{av}}^{}`$ in Eq. (8), and $`\mathrm{Var}(q_l)`$ in Eq. (10) tend to a non-zero constant in RSB, whereas they tend to zero with a power of $`L`$ in TNT and DP.

## III Results at Zero Temperature

We first discuss the results for the constrained average of $`1q_l`$, including only samples $`|q|0.5`$, since this yields $`dd_s`$ independent of $`\theta ^{}`$, see Eq. (8). The results are shown in Fig. 1. Note that the data only depend slightly on $`ϵ`$, indicating only small deviations from the expected scaling form.

The results of a DP/TNT fit to $`[1q_l]_{\mathrm{av}}^{}=L^{(dd_s)}(a+bL^c)`$ are presented in Table 3, while the corresponding RSB fits to $`[1q_l]_{\mathrm{av}}^{}=a+bL^c`$ are shown in Table 4. In the DP/TNT fits we find that $`dd_s`$ is close to zero.

Both DP/TNT and RSB fits are acceptable ($`\chi ^2/\text{ndf}1`$). However, in fits to a nonlinear model, one cannot convert $`\chi ^2/\text{ndf}`$ to a confidence limitPress et al. (1995) even if the data have a normal distribution. Similarly, unlike for the case for fits to a linear model, the error bars do not correspond to a $`68\%`$ confidence. We are particularly interested to get a confidence limit on the value of $`dd_s`$ in the DP/TNT fits. We do this by computing $`\chi ^2`$ as a function of $`dd_s`$, minimizing with respect to the other parameters ($`a`$, $`b`$, and $`c`$). The probability of the fit $`P`$ is proportional to $`\mathrm{exp}(\chi ^2/2)`$ which we numerically integrate to get the cumulative probability for $`x=dd_s`$:

$$Q(x)=^xP(x^{})𝑑x^{}.$$
(11)
The results are shown in Fig. 2. The data for $`ϵ=0.5`$ and $`2.0`$ constrain $`dd_s`$ to zero or a small positive value. The data for $`ϵ=1.0`$ constrain $`dd_s`$ less and allow a range of negative values which are unphysical. At a $`68\%`$ confidence level the data are consistent with

$$0dd_s0.05,$$
(12)
apart from the $`ϵ=0.5`$ data which would exclude zero at the $`68\%`$ level but, from Fig. 2, are consistent with it at the $`86\%`$ level. We take Eq. (12) to be our estimate for $`dd_s`$. It is consistent with the RSB prediction of zero and also consistent with a small non-zero value in the DP/TNT scenarios.

Data for the average of $`1q_l`$ over all samples are shown in Fig. 3 along with RSB fits. The data show curvature indicative that corrections to scaling have to be included. The DP/TNT fits, Eq. (6), are shown in Table 5, while the RSB fits (which fix $`\mu _l`$ to zero) are shown in Table 6. Both fits have acceptable $`\chi ^2`$.

The cumulative probabilities shown in Fig. 4 give a lot of weight to unphysical negative values of $`\mu _l`$, especially for $`ϵ=1`$. For all values of $`ϵ`$ the weight is small for $`\mu _l`$ greater than about $`0.10`$ so we conclude that

$$0\mu _l0.10.$$
(13)
We should, perhaps, be cautious about this statement in view of the large weight at negative values of $`\mu _l`$ in Fig. 4. However, Eq. (13) is consistent with Eq. (12) and the result of Ref. Katzgraber and Young, 2003 that $`\theta ^{}0`$.

## IV Results at Finite Temperature

In this section we study the model at temperatures well belowKatzgraber and Young (2003) $`T_c0.63`$. Figure 5 shows data for the variance of the link overlap for several low temperatures. The data show strong curvature indicative that a simple fit of the form $`aL^{\mu _l}`$ is improbable and corrections to scaling must be included.

Fits to Eq. (10) (DP/TNT picture) are shown in Table 7, and fits to the RSB picture (in which $`\mu _l`$ is fixed to be zero) are shown in Table 8. The quality of the fits is acceptable.

However, the DP/TNT fit for $`T=0.05`$ gives an unphysical negative value for $`\mu _l`$ with a very small amplitude $`a`$. To clarify this situation, we plot, in Fig. 6, $`\chi ^2`$ as a function of $`\mu _l`$, optimizing with respect to the other parameters ($`a,b`$, and $`c`$). For $`T=0.05`$, $`\chi ^2`$ is quite small out to very large negative values of $`\mu _l`$ (not shown) and increases rapidly for $`\mu _l`$ greater than about $`0.12`$. Since physically $`\mu _l`$ cannot be negative, the only conclusion we can deduce from the $`T=0.05`$ data is that $`\mu _l`$ lies between zero and about $`0.12`$, consistent with the result from the $`T=0`$ data in Eq. (13). The data for $`\chi ^2`$ for $`T=0.10`$ in Fig. 6 has a minimum at $`\mu _l=0.10`$ but it is shallow and $`\mu _l=0`$ has only a slightly greater $`\chi ^2`$ value. The $`T=0.10`$ data are therefore also consistent with Eq. (13). The data at higher temperatures, $`T=0.16`$ and $`0.23`$ have a deeper minimum at nonzero $`\chi ^2`$, suggesting that $`\mu _l=0`$ is somewhat unlikely, but experience from short-range systemsKatzgraber et al. (2001) suggests that estimates of $`\mu _l`$ at finite $`T`$ are effective exponents which need to be extrapolated to $`T=0`$ to get close to the asymptotic value. Hence we do not feel that the results at $`T=0.16`$ and $`0.23`$ rule out $`\mu _l=0`$.

Overall, the finite-$`T`$ data are consistent with $`\mu _l`$ in the range given by Eq. (13) which came from the $`T=0`$ data, and do not constrain $`\mu _l`$ any further.

## V Conclusions

We have studied the geometry of the large-scale, low-energy excitations in a one-dimensional Ising spin glass where the interactions fall off as $`r^\sigma `$ with $`\sigma =0.75`$, both at $`T=0`$ and at temperatures well below the spin-glass transition temperature. We find that the fractal dimension of the surface of these excitations, $`d_s`$, lies in the range $`0dd_s0.05`$. This is consistent with the RSB picture ($`dd_s=0`$). It is also consistent with the DP/TNT picture ($`dd_s>0`$) but with a small value of $`dd_s`$. Substantial corrections to scaling had to be incorporated into all the fits.

We have also estimated the exponent $`\mu _l=\theta ^{}+2(dd_s)`$, where $`\theta ^{}`$ characterizes the dependence of the energy of droplet excitations on their length scale. We find it to be in the range $`0\mu _l0.10`$, which is consistent with the value for $`dd_s`$ in Eq. (12) and our earlier resultKatzgraber and Young (2003) that $`\theta ^{}0`$. Note that this result for $`\theta ^{}`$, if also valid in the thermodynamic limit, is inconsistent with the DP.

By studying a one-dimensional model, we have been able to study a much larger range of sizes, $`16L512`$, than is generally possible in spin glasses. However, in the absence of a good understanding of corrections to scaling in spin glasses, we still cannot rule out the possibility that different behavior may occur in the thermodynamic limit.

###### Acknowledgements.

We would like to thank K. Tran for carefully reading the manuscript. The simulations were performed on the Asgard cluster at ETH Zürich. We are indebted to M. Troyer and G. Sigut for allowing us to use the idle time on the Asgard cluster. The work of APY is supported by the National Science Foundation under Grant No. DMR 0086287.


text:
# Precision Primordial 4He measurement from the CMB

## I Introduction

Big-bang nucleosynthesis (BBN) and the cosmic microwave background anisotropy (CMB) are two pillars of the hot big bang model.

The theory of BBN has long stood as an emblem of the predictive power of the Big Bang model wssok ; sarkar ; osw ; PDG\_reveiw . BBN has long provided the most reliable measurement of the cosmological baryon density. However, the CMB is rapidly becoming the preferred method for determining the baryon density, with its rapidly increasing precision. With the CMB, light element abundance observations become a powerful probe of the universe.

The observation and analysis of cosmic microwave background (CMB) anisotropies have attracted a great deal of attention in recent years due to their unique relevance for cosmological theory (see DodelsonHu  for a recent review). A flood of observational results have been published during the last two yearsToco98 ; DASI ; Boomerang02 ; Maxima ; VSA ; ACBAR ; CBI ; WMAP . These observations taken together measure CMB anisotropies over a large range of angular scales. The CMB is sensitive to the properties of the photon-baryon fluid and hence allows a precision determination of the baryon density at redshift $`z1000`$.

It is therefore apparent that combining BBN and CMB provides an opportunity for meaningful consistency checks on the standard cosmology and has the potential to be powerfully predictive probe of nuclear and particle astrophysics at low and high redshift Schramm\_Turner\_98 ; CFO2 ; CFO3 ; barger ; Hansen\_01 ; Kneller\_01 ; DiBari\_Foot\_01 . Given the constraint on the baryon density from the CMB, BBN yields a tight prediction of primordial $`{}_{}{}^{4}\mathrm{He}`$ abundance.

We also explore the promise of combining the CMB data with measurements of the deuterium (D/H) abundance, showing that current measurements of (D/H) can combine with the CMB constraints to reduce the error bars on $`{}_{}{}^{4}\mathrm{He}`$ by another factor of 2.

In section 2 we explain the data (CMB, BBN) and the method used to determine the likelihood surface in parameter space - the Metropolis-Hastings Markov Chain Monte Carlo (MCMC) algorithm. In section 3 the parameter confidence intervals that are extracted from that likelihood surface are discussed. In section 4 we discuss the implications of our results and how precision might be increased with further cosmological data.

## II methods

Our cosmological parameters are measured by use of the Markov Chain Metropolis Algorithm. We allow the 4He mass fraction to float as an independent variable, yielding the following parameter space: $`\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda },n_s,h,\mathrm{\Omega }_Bh^2,Y_p,\tau ,r,n_t`$ where $`\mathrm{\Omega }_{(m,\mathrm{\Lambda },B)}`$ is the matter, cosmological constant, and baryon contents, respectively, $`Y_p`$ is the fraction of the baryons in 4He by mass, $`h`$ is the Hubble constant in units of $`100km/sMpc`$, $`n_{(s,t)}`$ are the power-law index of the primordial scalar and tensor perturbations respectively, $`r`$ is the fraction of the observed CMB quadrupole that is tensor and $`\tau `$ is the optical depth to the last-scattering surface (that is, re-scattering of CMB photons by reionization is allowed for). Note that in our selection of the parameter space simplifying assumptions have been made: adiabatic, scale-free primordial perturbations, the universe contains only cold dark matter and a cosmological constant and the neutrino species are strictly those of the Standard Model. To determine the region of this parameter space allowed by experimental data, one must sample the space over a wide range of points. At a given point the relative likelihood of the parameter values yielding the observations must be determined, and the range of points sampled must adequately cover the space. Our primary likelihood calculation is a comparison of the simulated CMB spectra produced by CMBFAST CMBFAST  against the WMAP CMB experiment WMAP , along with smaller-scale (bin $`\mathrm{}_{eff}>350`$) data from the following experiments: Toco98 Toco98 , DASI DASI , Maxima Maxima , VSA VSA , ACBAR ACBAR , Boomerang02 Boomerang02 , and CBI CBI  111For WMAP, the published likelihood function was used. For the other experiments BJK formalism was used, with the BJK parameters obtained from either RADPACK, or the collaborations directly.. We include the published calibration uncertainties for each experiment and find the maximum likelihood value for these parameters at each point in the cosmological parameter space.

Figure 1 shows the data we used as well as the most likely 68% of our inferred theory power spectra.

One could attempt to sample the parameter space on a uniform grid, but the high dimensionality coupled with the computational demands of CMBFAST makes this impossible in a reasonable amount of time. Instead, we implemented the Metropolis-Hastings MCMC AlgorithmMetropolis ; Knox ; BridleLewis : starting from the current point in parameter space $`\stackrel{}{X}_i`$ one proposes a test point $`\stackrel{}{p}`$, drawn from a distribution described by a density function called the proposal density: $`P(\stackrel{}{p}|\stackrel{}{X}_i)`$. The choice of proposal density is somewhat arbitrary, but a poor choice will cause the parameter estimation procedure to be inefficient. The likelihood relative to $`\stackrel{}{X}_i`$ is computed at the test point $`\stackrel{}{p}`$. If the likelihood at $`\stackrel{}{p}`$ is greater, then point $`\stackrel{}{X}_{i+1}=\stackrel{}{p}`$. Otherwise there is a probability that $`\stackrel{}{X}_{i+1}=\stackrel{}{p}`$ equal to the likelihood of $`\stackrel{}{p}`$ divided by the likelihood of $`\stackrel{}{X}_i`$. If the proposal density is not symmetric in its arguments then this probability needs to be corrected by the factor $`P(\stackrel{}{X}_i|\stackrel{}{p})/P(\stackrel{}{p}|\stackrel{}{X}_i)`$ in order to enforce detailed balance. If $`\stackrel{}{X}_{i+1}`$ is not set to the the point $`\stackrel{}{p}`$, then it is set to the current point: $`\stackrel{}{X}_{i+1}=\stackrel{}{X}_i`$. After a sufficient number of iterations the resulting density of points $`\left\{\stackrel{}{X}_i\right|i=1\mathrm{}n\}`$ asymptotically approaches the likelihood function on parameter space. What constitutes a “sufficient” number of points is in general difficult to determine (and impossible to determine with absolute certainty) - though there are tests which are good indicators. We use a test for our Markov Chain which was suggested by Gelman92  and also used by the WMAP teamWMAP\_Verde .

We tune our Markov Chain code for optimal efficiency by first finding the maximal likelihood point in parameter space using a global maximization method (within our prior space, see below) then using it as the starting point for a sample Markov Chain. The variance of the sample chain is used to compute a step size matrix that will be used by the main chains. An efficient Markov Chain should take steps that are not too large or too small - either will result in an inefficient, slowly converging chain.

We choose our proposal density to be a multivariate Gaussian, whose covariance matrix $`V_s`$ is proportional to the covariance of the sample chain $`V_a`$. The proportionality constant $`\alpha (D)`$ is chosen so if the underlying distribution were Gaussian, $`50\%`$ of the Markov chain points would be accepted. Its value depends on the number of dimensions of parameter space, $`D`$ (for us $`D=9`$). We found:

$$\alpha \left(D\right)0.54784D0.36159$$
(1)
We ran two sets of Markov Chains, $`A`$ and $`B`$, each consisting of $`20`$ independent chains. Each chain was started at a point chosen from the distribution of the sample chain. Set $`A`$ had only very weak top-hat priors ($`\mathrm{\Omega }_Bh^2[0.014,0.030],`$ $`Y_p[0.13,0.34]`$, $`h[0.45,0.95]`$, $`\mathrm{\Omega }_m[0.03,1.00]`$, $`\mathrm{\Omega }_\mathrm{\Lambda }[0.00,0.97]`$, $`n_s[0.5,1.4]`$, $`n_t[3.0,3.0]`$, $`r[0.0,3.0]`$, $`\tau [0.0,1.0]`$), whereas set $`B`$ additionally had two strong priors: a BBN consistency condition between $`\mathrm{\Omega }_Bh^2`$ and $`Y_p`$ CFO1 , and the constraint $`h=0.74+0.110.094`$ Hproj  222The strong prior $`h=0.74+0.110.094`$ is an asymmetric Gaussian with plus and minus sigmas as given, and represents the constraint coming from the Hubble key project Hproj . The BBN consistency condition is simply this: for a given baryon density one expects BBN to produce a certain abundance of 4He, with some theoretical error (mostly driven by uncertainties in the nuclear cross-sections). In set $`A`$ we treat $`\mathrm{\Omega }_Bh^2`$ and $`Y_p`$ as two independent variables. In set $`B`$ we enforce theoretical self-consistency between those variables. This non-CMB constraint is incorporated in a Bayesian way. The additional information is included as an additional prior and the density to be sampled from (the posterior density) becomes

$$\begin{array}{c}P_{tot}\left(\stackrel{}{X}\right)=P_{CMB}\left(\stackrel{}{X}\right)P_{BBN}(\mathrm{\Omega }_Bh^2,Y_p)P_{Hubble}\left(h\right)\\ =P\left(\stackrel{}{X}|CMB\right)P\left(Y_p|\mathrm{\Omega }_Bh^2\right)P\left(h|Hubbledata\right)\end{array}$$
(2)
where $`P\left(x|y\right)`$ is the conditional density for getting $`x`$, given $`y`$. Note that $`P\left(Y_p|\mathrm{\Omega }_Bh^2\right)`$ is purely a theoretical prior enforcing the BBN relation, and contains no abundance measurements. The Markov Chain then automatically explores the new, more constrained region of parameter space. As one might expect, sets $`A`$ and $`B`$ differ significantly in their parameter space coverage, and thus their proposal densities and chain starting points were determined separately.

Though a Markov Chain approach saves significant computational time, it is difficult to guarantee after some number of points that the chain has converged sufficiently to the true, underlying distribution. Indeed, a chain cannot tell one anything about a region it has not yet visited. We use a convergence test suggested by Gelman & Rubin Gelman92 , which was also employed by WMAP WMAP\_Verde .

We have generalized this criterion to multiple dimensions, keeping in mind that any convergence test must be covariant (if a transformation of parameter space can change the determination of “convergence”, then the test is a bad one). Each chain out of the $`20`$ has its own mean and variance 333For $`D>1`$ we use the word variance to mean covariance matrix.. If each chain reflected the underlying distribution, then the variance of the means of the chains should be much less than the variance of the underlying distribution. We thus compute the variance of the chain means, multiply with the inverse of the average chain variance, and take the trace:

$$\begin{array}{c}U\frac{1}{N1}\underset{j=1}{\overset{N}{}}\left(\stackrel{}{\overline{X}}_j\stackrel{}{\overline{X}}\right)\left(\stackrel{}{\overline{X}}_j\stackrel{}{\overline{X}}\right)\\ W\frac{1}{N}\underset{j=1}{\overset{N}{}}\frac{1}{n_j1}\underset{i=1}{\overset{n_j}{}}\left(\stackrel{}{X}_{j,i}\stackrel{}{\overline{X}}_j\right)\left(\stackrel{}{X}_{j,i}\stackrel{}{\overline{X}}_j\right)\\ \mu Tr\left[UW^1\right]/D\end{array}$$
(3)
where $`N`$ is the number of chains (20), $`n_j`$ is the number of points in chain $`j`$, $`\stackrel{}{\overline{X}}`$ is the total mean, $`\stackrel{}{\overline{X}}_j`$ is the mean of chain $`j`$. We require that $`\mu <0.1`$. Set $`B`$ easily satisfies this criteria with $`30,000`$ points, whereas set $`A`$ required about $`60,000`$ points. The average chain variance, $`W`$, is used because this underestimates the variance of the distribution until convergence is attained.

As a self-consistency check, one can take the point distribution of set $`A`$ and combine it with the BBN $`\mathrm{\Omega }_Bh^2Y_p`$ and Hubble Key project priors. The resulting distribution should be the same as set $`B`$. The extent to which these distributions differ is a measure of non-convergence of the sets. We determined that the $`68\%`$ confidence regions of these distributions more than $`95\%`$ overlap in the $`(\mathrm{\Omega }_Bh^2,Y_p)`$ plane. A Markov Chain can be combined with a prior after the generation of the chain by assigning a weight to each point. The likelihood of a region in parameter space is then the weighted density of the points in that region. Because a Markov Chain maintains the full $`D`$-dimensional likelihood distribution in the parameter space, after it is generated the chain may be convolved with any arbitrary other likelihood function in that parameter space. Thus one can generate a Markov Chain distribution for WMAP alone, and chose any subset of the other cosmological datasets to convolve it with - for very little additional CPU cost. This is the basis of the *Cosmic Concordance Project* (CCP) web-site http://galadriel.astro.uiuc.edu/ccp/ where the parameter constraints from chains $`A`$ and $`B`$ can be explored and combined with other cosmological datasets and priors. Further details about the CCP and about our parameter estimation methodology will be given in CCP .

## III Results

As our Markov Chain sets $`A`$ and $`B`$ approached $`60,000`$ and $`30,000`$ points respectively, convergence test eq. 3 gave $`\mu 0.050.06`$ 444For this value we concatenated $`5`$ chains into $`1`$, thus transforming $`20`$ chains per set to $`4`$. and we declared our chains sufficiently converged to provide reliable statistics. Figure 2a shows the 68% confidence region of set $`A`$ in the $`(\mathrm{\Omega }_Bh^2,Y_p)`$ plane. Thus from CMB data alone we find $`\mathrm{\Omega }_Bh^2=0.025_{0.0026}^{+0.0019}`$ and $`Y_p=0.250_{0.014}^{+0.010}`$

It is worth noting also that Markov chain set A yields a tight constraint on the neutrino number during BBN: Allowing the number of BBN neutrinos to float, one can find a 68% confidence interval for the number of BBN neutrinos needed to yield the determined $`Y_p`$ from the determined $`\mathrm{\Omega }_Bh^2`$. Given a BBN code that computes the probability density $`P(Y_p|\mathrm{\Omega }_Bh^2,N_{\nu ,eff})`$ (where the stochasticity is due to the measurement uncertainty in the relevant nuclear cross-sections) we can compute $`P(N_{\nu ,eff})`$ based on CMB data and BBN theory as

$`P(N_{\nu ,eff}|\mathrm{BBN},\mathrm{CMB})`$ $`=`$
$`{\displaystyle 𝑑Y_pd(\mathrm{\Omega }_Bh^2)P(N_{\nu ,eff}|Y_p,\mathrm{\Omega }_Bh^2)P(Y_p,\mathrm{\Omega }_Bh^2)}`$ $`=`$
$`{\displaystyle 𝑑Y_pd(\mathrm{\Omega }_Bh^2)P(Y_p|\mathrm{\Omega }_Bh^2,N_{\nu ,eff})P(Y_p,\mathrm{\Omega }_Bh^2)}.`$
The first term under each integral enforces the BBN relation and the second term the CMB posterior. The second equality holds true if we assume flat priors $`P(Y_p|\mathrm{\Omega }_Bh^2,BBN)=const.`$ and $`P(N_{\nu ,eff}|\mathrm{\Omega }_Bh^2,BBN)=const`$ for the range of parameter space of interest 555 For simplicity, we approximate $`P(Y_p,\mathrm{\Omega }_Bh^2|CMB)=P(Y_p|CMB)P(\mathrm{\Omega }_Bh^2|CMB)`$ for this calculation. We also assume that the CMB does not constrain $`N_{\nu ,eff}`$ (instead, we treat the effective number of neutrino species that are relevant for the CMB anisotropy as independent of $`N_{\nu ,eff}`$) and that there is no a priori preference for any value of $`N_{\nu ,eff}`$ or $`Y_p`$ for any value of $`\mathrm{\Omega }_Bh^2`$. . We find $`N_{\nu ,eff}=3.02_{0.79}^{+0.85}`$. $`N_{\nu ,eff}`$ is consistent with the standard model value of 3 and previous studies CFO2 ; CFO3 ; barger ; Hansen\_01 ; Kneller\_01 ; DiBari\_Foot\_01 .

However, $`\mathrm{\Omega }_Bh^2`$ and $`Y_P`$ are jointly constrained by BBN theory, and thus are not really independent variables. Adopting the standard BBN model ($`N_{\nu ,eff}(\mathrm{BBN})`$=3) of CFO1  yields a consistency relation between $`\mathrm{\Omega }_Bh^2`$ and $`Y_P`$. The $`68\%`$ confidence region of this consistency relation appears in Figure 2a as a narrow band (narrow enough that the upper and lower bounding curves appear to merge). Enforcing this condition greatly increases the precision of parameter estimation, as evident in table 1, with dramatic affect on $`Y_P`$ measurement: $`Y_P=0.2493_{0.0010}^{+0.0007}`$. This is simply a result of the accuracy with which $`{}_{}{}^{4}\mathrm{He}`$ is determined by BBN ($`0.1\%`$). In Figure 2b we have zoomed in on this CMB-BBN concordance region. Also shown is the $`68\%`$ confidence region of set $`B`$ shaded as solid red. As one would expect, set $`B`$ agrees with the product of the CMB (set $`A`$) and BBN (consistency band) likelihoods. It is important to note that in Figure 2a the agreement between the CMB and BBN allowed regions need not have happened. Instead, the BBN consistency band might not have passed through the high CMB likelihood region, which would have forced one to consider a BBN scenario other than the standard model one with 3 neutrinos. The CMB-BBN agreement reaffirms the standard BBN scenario.

Model A and Model B compare quite well to $`{}_{}{}^{4}\mathrm{He}`$ observations. Olive, Skillman and Steigman (1997) Olive\_Skillman\_Steigman\_97  and Fields and Olive (1998) Fields\_Olive\_98  find $`Y_p=0.238\pm 0.002`$, while Izotov and Thuan (1998) Izotov\_Thuan\_98  find $`Y_p=0.244\pm 0.002`$. The errors cited are statistical only. Comparing these numbers, not only are they discrepant from each other, but they lie below the mean value determined in this evaluation. However, Olive and Skillman (2001) Olive\_Skillman\_16  critically evaluate the methods used in determining $`Y_p`$ and find a lower bound to a systematic error, $`\sigma _{sys}0.005`$. This systematic error is added in quadrature with the statistical error to determine the total error, increasing the errors to $`0.0054`$. The $`Y_p`$ observations are brought into marginal accord with each other and the CMB; both lie systematically lower than the CMB determined value. As discussed earlier, the systematic error used is only a lower bound, and as such the true errors are most likely larger than those quoted.

In Figure 2 the allowed $`\mathrm{\Omega }_Bh^2`$ is very large. Any other data that can reduce the allowed $`\mathrm{\Omega }_Bh^2`$ range will have the additional benefit of refining the precision of the $`Y_P`$ measurement. As an example, we consider $`D`$ abundance in Figure 3. The value of $`D/H`$ is still a somewhat open question due to low number statistics, and thus we demonstrate the effects of two different $`D`$ abundances. In Figure 3a we use the average of the 2 multiple absorption line systems of OMeara01\_Dma1 ; Kirkman03\_Dma2 : $`D/H=\left(2.49\pm 0.18\right)\times 10^5`$. The blue shaded band is the $`68\%`$ confidence baryon density range allowed by this $`D/H`$ value as determined with the BBN theory of CFO1 . The solid red region is the $`68\%`$ confidence region resulted from the convolution of the $`D`$ baryon range and the concordance region of Figure 2b. Thus we find $`Y_p=0.2491_{0.0005}^{+0.0004}`$, $`\mathrm{\Omega }_Bh^2=0.0237_{0.0012}^{+0.0010}`$. Alternatively, one can use a conservative combination of 5 $`D/H`$ measurements, including the two multiple absorption line systems employed above OMeara01\_Dma1 ; Kirkman03\_Dma2 ; Burles98a\_Dwa1 ; Burles98b\_Dwa2 ; Pettini01\_Dwa3 : $`D/H=\left(2.78\pm 0.29\right)\times 10^5`$ (the overall error increase is because the other three systems are not consistent with each other or the multiple absorption line systems, a hint of an underlying systematic error for the single absorption line systems). For this $`D`$ abundance we find $`Y_p=0.2488_{0.0005}^{+0.0004}`$ and $`\mathrm{\Omega }_Bh^2=0.0230_{0.0012}^{+0.0008}`$.

## IV Conclusions

This work has been based on two general ideas; (1) BBN and the CMB independently probe two different epochs, providing valuable consistency checks for the underpinnings of the Standard Cosmological Model; (2) having established that the cosmological model agrees remarkably well with these very different observational probes we use these data to to make a precision measurement of the 4He abundance. We have presented an analysis of all recent CMB data, in which we have determined the cosmic baryon density and the primordial helium abundance. We found $`\mathrm{\Omega }_Bh^2=0.0250_{0.0026}^{+0.0019}`$ and $`Y_p=0.250_{0.014}^{+0.010}`$ at $`68\%`$ from CMB data alone. This is consistent with 3 standard model neutrinos during BBN.

We have shown that this is fully consistent with the predicted 4He abundance from BBN theory, and marginally consistent with 4He observations. The likely source of this slight discrepancy is an underestimate of the dominant, systematic uncertainties in the 4He observations, which now seems affirmed with the CMB determination of $`Y_p`$. The agreement between the CMB-only set $`A`$ confidence region of $`\mathrm{\Omega }_Bh^2`$, $`Y_p`$ and the consistency band based on BBN theory shown in Figure 2a reaffirms the standard BBN model. Thus using BBN theory, we can effectively remove $`Y_p`$ as a freely floating variable, enforcing the $`\mathrm{\Omega }_Bh^2`$-$`Y_p`$ BBN relation in CMB data analysis. Given this, we found the incorporation of BBN theory into parameter extraction from CMB data results in a precision measurement of the $`{}_{}{}^{4}\mathrm{He}`$ abundance. We find $`\mathrm{\Omega }_Bh^2=0.0245_{0.0029}^{+0.0015}`$ and $`Y_p=0.2493_{0.0010}^{+0.0007}`$.

Using this CMB-BBN determined $`Y_p`$, one can study and constrain stellar evolution Bono\_02 ; Cassisi\_Salaris\_Irwin\_03 . One can similarly study the nucleosynthetic history of all the light element abundances as discussed in CFO3  and references therein.

We show the promise of incorporating deuterium abundance observations, yielding $`Y_p=0.2491_{0.0005}^{+0.0004}`$, $`\mathrm{\Omega }_Bh^2=0.0237_{0.0012}^{+0.0010}`$ or $`Y_p=0.2488_{0.0005}^{+0.0004}`$, $`\mathrm{\Omega }_Bh^2=0.0230_{0.0012}^{+0.0008}`$ depending on which systems are used to measure the deuterium abundance.

The addition of the $`D`$ abundance observations are only one example of many possible cosmological datasets that might be incorporated into parameter extraction to increase precision. An experiment may constrain a parameter directly, or may reduce degeneracy in a related parameter. For example, using large-scale structure information to reduce the residual $`n_s`$-$`\mathrm{\Omega }_Bh^2`$ degeneracy in the current CMB data would also increase the precision of the $`Y_P`$ determination. Also, further light element observations and CMB anisotropy experiments will refine this concordance and sharpen BBN and the CMB as tools for precision cosmology. Due to the effect of the $`{}_{}{}^{4}\mathrm{He}`$ abundance on the damping tail, this may improve the constraint on a possible running of scalar spectral index. These are the topics of on-going work and can be further explored at the *Cosmic Concordance Project* web-site: http://galadriel.astro.uiuc.edu/ccp/ CCP .

Many thanks to Prof. Brian Fields for helpful discussions. The work of R.H.C was provided by the National Science Foundation grant AST-0092939. This work was partially supported by National Computational Science Alliance under grant AST020003N and utilized the NCSA Platinum system.