text:
# Evidence for Planet-induced Chromospheric Activity on HD 179949

## 1 Introduction

Current planet detection methods provide basic information: a minimum mass if the orbital inclination is not known and an albedo-dependent estimate of surface temperature. In the case of the transiting system, HD 209458, the planet’s density is also determined with a possible atmospheric sodium detection (Charbonneau 2002) and excess Lyman $`\alpha `$ absorption during transit (Vidal-Madjar et al. 2003). There still remains a lack of constraints on extrasolar planetary structure leaving astronomers to explore new observational probes.

Cuntz, Saar & Musielak (2000) suggested that there may be an observable interaction between a parent star and a close-in giant planet, specifically an external heating of the star’s outer atmosphere. The effect could be tidal, magnetic, or a combination of the two. In the case of tidal heating, the acoustic wave energy thought to contribute to the heating of the upper atmosphere has a dependence on the local turbulent velocity ($`v_t`$) as $`v_t^8`$ (Muzielak et al. 1994). The steepening of the density gradient in the two tidal bulges would increase the turbulent velocity producing excess shocks, waves, and flows in the upper atmosphere of the star. Magnetic wave energy has a $`v_t^6`$ dependence (Ulmschneider & Muzielak 1998) such that an increase in $`v_t`$ would enhance local dynamo generation and hence the local magnetic energy will also increase. Thus, even small increases in $`v_t`$ would cause a dramatic intensification in nonradiative heating. If planet-induced heating is confined to a narrow range in stellar longitude, the heated regions would track the planet. This implies that the period of any observed activity would be correlated with the planet’s orbit such that tidally induced activity has a period of $`P_{orb}/2`$ and magnetic activity, a period of $`P_{orb}`$. It is also possible that the increased dynamo action would be confined to the turbulent layer below the convection zone (the tachocline). In this case, any dynamo action would be spread out over a wide range of longitudes and the excess heating would not show strong periodic variation (Saar & Cuntz 2001).

Since the tidal and magnetic interaction depend on the distance from the planet to the star as $`1/d^3`$ and $`1/d^2`$, respectively, it is best to look at the tightest systems. Of the $``$ 100 known extrasolar planets, 16% have semi-major axes of less than 0.1 AU and masses comparable to Jupiter’s (Schneider 2003). It is expected that these planets have magnetic fields similar to Jupiter’s (4.3 G). However if tidally locked, the planets’ fields might be substantially smaller. It is also reasonable to assume that the magnetic interaction would be greatest in the outermost layers of the star, namely the chromosphere, transition region and the corona due to their proximity to the planets and their nonradiative heat sources. Our program stars have orbital periods between 3.1 and 4.6 days, eccentricities $``$ 0 and semi-major axes $`<0.06`$ AU. These systems offer the best chance of observing upper atmospheric heating. The five systems we chose were $`\tau `$ Boo, HD 179949, HD 209458, 51 Peg and $`\upsilon `$ And and their planetary system parameters are listed in Table 1.

## 2 Chromospheres and Ca II H & K Emission

Temperatures in the chromosphere can reach 20,000 K over many scale heights implying that the temperature gradients are not very high. This is due to the large opacity of resonance lines such as H I, Ca II, Mg II and the Lyman continuum which then in the low-density gas radiate and cool the chromosphere. Most chromospheric indicators are UV emission lines not easily accessible from ground-based telescopes. However, the optical resonance lines of Ca II H & K at 3933 Å and 3968 Å exhibit chromospheric emission in late-type stars (later than F0) and are sensitive probes of temperature and electron density. An example of HD 179949 Ca II H & K reversals are shown in Figure 1. The double-peaked emission is formed in the chromosphere at T $``$ 8000 K where the source function is essentially monochromatic. The central emission is produced at a temperature of $``$ 20,000 K (Montes et al. 1994) which appears as a self-reversal due to the large optical depth of the line. The asymmetry of the emission in late-type dwarfs is thought to be a consequence of upflows in the chromosphere’s plage-like regions. (See Linsky 1980 for a comprehensive observational and theoretical review.)

The broad, deep photospheric absorption and the low level background continuum allow the reversal to be seen at higher contrast. Because of these and the accessibility of the lines from the ground, the Ca II H & K reversals are an optimal choice with which to monitor chromospheric heating for the sun-like stars known to host short-period giants planets. Saar and Cuntz (2001) looked at the Ca II IR line at 8662 Å, an interlocking line with the Ca II H & K resonance lines, with spectra of resolution R = 50,000 and S/N typically of 200. They looked at five stars with close-in giant planets (including $`\tau `$ Boo, $`\upsilon `$ And, and 51 Peg) by simulating the IR equivalent to the Mt. Wilson $`S_{HK}`$ index and found nothing at the 3 to 5% level.

## 3 Observations

High-resolution (R = 110,000) and high signal-to-noise (S/N $``$ 500 per pixel in the continuum and 150 in the H & K cores) spectra were obtained on 3 observing runs at the Canada-France-Hawaii Telescope (CFHT) with the fiber-fed coudé spectrograph, $`Gecko`$. The single spectrum spanned 60 Å and was centered on 3947 Å. A representative spectrum of HD 179949 is shown in Figure 1. Five nights were allocated in each of three semesters: August 2001, July 2002 and August 2002. Poor weather and technical issues limited usable data to $``$ 10 nights. Spectra with comparable S/N were taken of two stars known not to have close-in giant planets, $`\tau `$ Ceti and the Sun (sky spectra were taken at dusk).

The data reduction consisted of the subtraction of the appropriate darks from the stellar, arc and flat-field exposures rather than using the convential biases to remove the baseline from each observation. The flats were combined and normalized to a mean value of 1 along each row of the dispersion axis. All the stellar exposures for each night were averaged in order to define a single aperture to use for the spectral extraction. This included a subtraction of the residual background level between orders. The aperture was used to extract 1-dimensional spectra from the stellar, comparison and flat-field exposures. A mean, normalized 1-dimensional flat was divided into the extracted spectra to obtain the most consistent flat-fielded spectra possible. The wavelength calibration was done with Thorium-Argon arcs. Heliocentric and differential radial velocity corrections were applied.

In 2001 Th-Ar arcs were taken with every target while in the 2002 observing runs, each stellar exposure was bracketed with an arc. This allowed differential radial velocities ($`\mathrm{\Delta }`$$`RV`$) to be measured to better than 20 m/s. Detailed discussion on the $`RV`$ can be found in Walker et al. (2003). Current orbital ephemerides and hence accurate phases ($`\pm 0.03`$) were determined for each observation by fitting a least-squares sine function with the known radial velocity amplitude and orbital period. A pure sinusoid was sufficient since all five systems have eccentricity $`e0.05`$ (Schneider 2003). Table 2 contains the July 2002 times of sub-planetary position along with the revised periods.

For each star, portions of the spectra, each 7 Å wide, were extracted centered on Ca II H, K and the strong Al I line at 3944 Å. As a means of normalizing each sub-spectrum, the continuum was set to 1 at the ends and fit with a straight line. The spectra were grouped by date and a nightly mean was computed for each of the three lines. An overall average was then taken over all the nights from which we measured nightly residuals. Each residual spectrum also had a broad, low-order curvature removed.

## 4 Results

The program stars were originally chosen for Doppler planet searches because of their weak chromospheric activity as indicated by the modest Ca II H & K emission flux measured by the Mt. Wilson $`S_{HK}`$ index. However there are marked differences in the reversal structure of each of these targets caused by their different rotational rates (which produce a range of emission strengths and Doppler broadening) as well as by other intrinsic mechanisms, presumably acoustic waves which depend on $`T_{eff}`$ and $`g`$. The integrated flux of the mean, normalized Ca II K core, including the photospheric contribution and the chromospheric emission, is listed in Table 1 and plotted against $`v`$sin$`i`$ in Figure 2. Even with the small number of program stars, there is a reasonable linear relationship between the two parameters as reported in Pasquini et al. (2000). HD 179949 lies significantly above the best-fit line. The Sun does as well, but this is known to be an effect of increased activity observed in the 2002 observing runs.111We observed the naked-eye sunspot grouping of July 2002 in transit during our CFHT run with the Sunspotter® Solar Telescope supplied by the Visitor Information Center at the Onizuka Center for International Astronomy.

The residuals of the normalized spectra (smoothed by 21 pixels) taken from the mean were used to generate the Mean Absolute Deviation (MAD = $`N^1\mathrm{\Sigma }|data_imean|`$ for $`N`$ spectra) plot shown in Figure 3 (top) for HD 179949. Superimposed on the plot is the corresponding K-core. The nightly residuals used to generate the MAD plot for HD 179949 are displayed just below. The activity seen in four of the five stars with planets, HD 179949, $`\tau `$ Boo, $`\upsilon `$ And and HD 209458, show significant deviations in the H & K reversals. The widths of the MAD plots (Figure 4) are the same as those of the reversals themselves as illustrated in Figure 3 for HD 179949. In the Sun’s case, images from the Michelson Doppler Imager (MDI) aboard the Solar Heliospheric Observatory (SOHO)222http://sohowww.nascom.nasa.gov/ and http://soi.stanford.edu/ confirm that the sunspot group was receding which caused the apparent asymmetry of the solar MAD plot. The strong coherence between the H and K activity can be seen in Figure 5 where the integrated mean absolute deviation for H and K is plotted for each star. The strong aluminum line shows no comparable variations. For all the program stars, the mean absolute deviation of the Al I line at 3944Å is less than 0.0005 of the normalized mean. This indicates that the photospheres of the stars are stable to this level providing additional evidence that the variations seen in the Ca II H & K reversals are not caused by upper-level photospheric changes but rather are confined to the chromosphere.

## 5 HD 179949: A Convincing Case of Planetary-induced Modulation

HD 179949 shows the greatest nightly modulation of all the stars in all three observing runs. The integrated flux of the residuals shown in Figure 3 (bottom) are plotted against orbital phase in Figure 6. The emission clearly increases by $``$ 4% when the planet is in front of the star, at phases near the sub-planetary point ($`\varphi 0`$) and is least when the planet is behind the star ($`\varphi `$ 0.5). The amplitude and period of this activity has persisted during the year between observations (equal to 108 orbits or at least 37 stellar rotations). A periodogram of even these few data points (Stellingwerf 1978) produces the most significant peak with a period of 3 days.

### 5.1 Rotation of HD 179949

HD 179949 has no photometrically measured rotation period so, as derived from the $`v`$sin$`i`$, $`P_{rot}`$ 9 days. The solid line in Figure 6 is a truncated, best-fit sine curve with $`P=P_{orb}`$ corresponding to the change in projected area of a bright spot on the stellar surface before being occulted by the stellar limb. The peak of the model leads the sub-planetary point by 0.17 in phase. If this effect is real, it could provide a constraint on future models of the interaction. Since Ca II emission is optically thin outside the self-absorption core, limb brightening needs to be included in the model. According to solar observations discussed by Engvold (1966), limb brightening would increase emission by 3%, essentially not changing the geometric model.

The best-fit spot model requires the bright spot to be at a latitude of $`30^{}`$ and a stellar inclination angle $`i`$ of $`87^{}`$. Since there were no transits detected for this system, an upper limit for the orbital inclination is set at $`83^{}`$ (Tinney et al. 2001). The dashed model in Figure 6 has $`i=83^{}`$ as it is generally assumed that the orbital and stellar inclination angles are comparable. These models are inconsistent with the possibility that the star may be tidally synchronized with the planet’s orbit ($`P_{rot}=P_{orb}`$). However, from Figure 2, the Ca II K emission in HD 179949 is unusually strong compared to the other targets. This implies that either the inclination of the star is low, $`21^{}`$ ($`v`$ 17.5 km s-1 and $`P_{rot}P_{orb}`$ 3 d.) or that $`P_{rot}`$ is closer to 9 days and the Ca II emission is significantly enhanced in HD 179949 relative to the other stars.

If the star is indeed tidally synchonized with the planet, as is believed to be the case for $`\tau `$ Boo, then the modulation in the emission observed with the 3.093-day period may be tracking a persistent stellar surface feature. However, tidal synchronization of HD 179949 is not a favoured scenario for two reasons. Firstly, the tidal synchonization timescale for the system $``$ 700 Gyr using a $`M_p`$sin$`i=0.84M_J`$. (See Equation 1 of Drake et al. 1998). Even if $`i=21^{}`$, as would be the case for $`v`$ = 17.5 km s-1 and $`M_p=2.34M_J`$, then the sychronization timescale is still very long at 79 Gyr. Secondly, ROSAT X-ray data list HD 179949 as having at least double the X-ray flux (a measurement independent of $`i`$) as compared to other single F8-9 dwarfs and $``$ 10 times that of the Sun’s. Hünsch et al. (1998) report the X-ray flux for HD 179949 to be $`L_X`$ = 41.0 $`\times 10^{27}`$ erg/s while for 36 UMa, a nearly identical star, $`L_X`$ = 18.4 $`\times 10^{27}`$ erg/s .

## 6 Summary

From the sample of stars observed, those with planets (with the exception of 51 Peg) show significant nightly variation in their Ca II H & K cores. $`\tau `$ Ceti which has no close planet remained very steady over all three observing runs. The Sun’s Ca II emission was stable within each run. However, a large naked-eye sunspot grouping appeared in the summer of 2002 and was detected as an enhancement in the H & K reversals accounting for its MAD signal (Figure 4). There is enhanced nightly activity in the stars with close-in giant planets whose variations may be cyclical and phase shifted. Further observations of these systems at varying orbital phases are required.

In the one clear case, HD 179949 exhibits repeated phase-dependent activity with enhanced emission near the sub-planetary point and less emission half an orbit later. This is consistent with a magnetic heating scenario and may be a first indirect glimpse at a magnetosphere of an extrasolar planet. Additional Ca II observations are crucial to confirm the continuity of the magnetic interaction as well as to establish better phase coverage. Due to the nearly exact 3-day orbit, it is impossible to get many distinct phases during one observing season from the ground. Space based observations will not only increase the phase coverage, but will also allow observations of other activity diagnostics such as transition region and coronal lines in the UV and FUV. Such mapping of planet-induced activity at varying stellar atmospheric depths will allow characterization and quantification of the physical interaction between the magnetically heated layers of the star and the planet’s magnetosphere.

Research funding from the Canadian Natural Sciences and Engineering Research Council (G.A.H.W. & E.S.) and the National Research Council of Canada (D.A.B.) is gratefully acknowledged. We are also indebted to the CFHT staff for their care in setting up the CAFE fiber feed and the Gecko spectrograph.


format: text

text:
# Determining Neutrino Mass from the CMB Alone

## Abstract

Distortions of Cosmic Microwave Background (CMB) temperature and polarization maps caused by gravitational lensing, observable with high angular resolution and high sensitivity, can be used to measure the neutrino mass. Assuming two massless species and one with mass $`m_\nu `$ we forecast $`\sigma (m_\nu )=0.15`$ eV from the Planck satellite and $`\sigma (m_\nu )=0.04`$ eV from observations with twice the angular resolution and $`20`$ times the sensitivity. A detection is likely at this higher sensitivity since the observation of atmospheric neutrino oscillations require $`\mathrm{\Delta }m_\nu ^2\begin{array}{c}>\hfill \\ \hfill \end{array}(0.04\mathrm{eV})^2`$.

Introduction. Results from the WMAP bennet03a  show the standard cosmological model passing a highly stringent test. With this spectacular success of the CMB as a clean and powerful cosmological probe, and of the standard model as a phenomenological description of nature, it is timely to ask what can be done with yet higher resolution and higher sensitivity such as offered by the Planck instruments and beyond. In this Letter we mostly focus on neutrino mass determination, with brief discussion of other applications.

Eisenstein et al. eisenstein99  found that the Planck satellite can measure neutrino mass with an error of 0.26 eV. This sensitivity limit is related to the temperature at which the plasma recombines and the photons last scatter off of the free electrons, $`T_{\mathrm{dec}}0.3`$ eV. Neutrinos with $`m_\nu \begin{array}{c}<\hfill \\ \hfill \end{array}T_{\mathrm{dec}}`$ do not leave any imprint on the last-scattering surface that would distinguish them from $`m_\nu =0`$.

Neutrinos with mass $`m_\nu \begin{array}{c}<\hfill \\ \hfill \end{array}T_{\mathrm{dec}}`$ would affect the amplitudes of gravitational potential peaks and valleys at intermediate redshifts. Massive neutrinos can collapse into potential wells when they become non-relativistic, while massless ones freely stream out. The observed galaxy power spectrum (which is proportional to the potential power spectrum at sufficiently large scales), combined with CMB observations can be used to put constraints on $`m_\nu `$ hu98a . At present such an analysis yields an upper bound on $`m_\nu `$ of $`0.3`$ eV spergel03 111No significant improvement is expected from combining Planck and the SDSS galaxy power spectrum hu98 ; eisenstein99 .

The alteration of the gravitational potentials at late times changes the gravitational lensing of CMB photons as they traverse these potentials seljak96 ; bernardeau97 . Including the gravitational lensing effect, we find that the Planck error forecast improves to 0.15 eV. We also show that more ambitious CMB experiments can reduce this error to $`0.04`$ eV. These mass ranges are interesting because the atmospheric neutrino oscillations require that at least one of the active neutrinos have $`m_\nu >0.04`$ to 0.1 eV. More detailed considerations beacom02  show that the sum of the active neutrino masses (which is what the CMB is most sensitive to) should be at least 0.06 eV.

Tomographic observations of the galaxy shear due to gravitational lensing can achieve sensitivities to $`m_\nu `$ similar to what we find here hu99 ; abazajian02b . Our work is distinguished by its sole reliance on CMB temperature and polarization maps which have different potential sources of systematic error. Complementary techniques are valuable since both of these will be very challenging measurements.

The most stringent laboratory upper bound on neutrino mass comes from tritium beta decay end-point experiments tritium  which limit the electron neutrino mass to $`\begin{array}{c}<\hfill \\ \hfill \end{array}2`$ eV. Proposed exeriments plan to reduce this limit by one to two orders of magnitude by searching for neutrinoless double beta decay ($`\beta \beta 0\nu `$) zdesenko03 . A Dirac mass would elude this search, but theoretical prejudice favors (and the see-saw mechanism requires) Majorana masses. Like the CMB and galaxy shear observations, these future $`\beta \beta 0\nu `$ experiments will be extremely challenging.

Lensing of the CMB. The intensity and linear polarization of the CMB are completely specified by the Stokes parameters, $`I`$, $`Q`$ and $`U`$ which are related to the unlensed Stokes parameters (denoted with a tilde) by $`X(𝐧)=\stackrel{~}{X}(𝐧+\delta 𝐧)`$ where $`X`$ stands for $`I`$, $`Q`$ or $`U`$. The deflection angle, $`\delta 𝐧`$, is the tangential gradient of the projected gravitational potential,

$$\varphi (𝐧)=2𝑑r\mathrm{\Psi }(r\widehat{𝐧},r)(rr_s)/(rr_s),$$
(1)
where $`r`$ is the coordinate distance along our past light cone, $`s`$ denotes the CMB last–scattering surface, $`\widehat{𝐧}`$ is the unit vector in the $`𝐧`$ direction and $`\mathrm{\Psi }`$ is the three-dimensional gravitational potential.

The statistical properties of the $`I`$, $`Q`$ and $`U`$ maps are most simply described in the transform space: $`a_T(𝐥)`$, $`a_E(𝐥)`$, and $`a_B(𝐥)`$ where $`a_T`$ is the spherical harmonic transform of $`I`$ and $`a_E`$ and $`a_B`$ are the curl–free and gradient–free decompositions, respectively, of $`Q`$ and $`U`$ kamionkowski97 ; seljak97 . In this transform space the effect of lensing by mode $`\varphi (𝐋)`$ (harmonic transform of $`\varphi (𝐧)`$) is to shift power from, e.g., $`\stackrel{~}{a}_T(𝐋𝐥)`$ to $`a_T(𝐥)`$. Lensing also mixes $`\stackrel{~}{a}_E`$ into $`a_B`$ and any $`\stackrel{~}{a}_B`$ into $`a_E`$ zaldarriaga98b , thus generating scalar B (curl) mode correlations.

Lensing smoothes out the features in the two-point functions, also called angular power spectra, $`𝒞_l^{\alpha \alpha ^{}}`$, where $`a_\alpha (𝐥)a_\alpha ^{}^{}(𝐥^{})=𝒞_l^{\alpha \alpha ^{}}\delta (𝐥𝐥^{})/[2\pi l(l+1)]`$ and $`\alpha `$ stands for $`T`$, $`E`$, or $`B`$ seljak96 . As explained later, in our analysis we use the unlensed power spectra, $`\stackrel{~}{𝒞}_l^{\alpha \alpha ^{}}`$. The information from lensing is added through the two-point function of the lensing potential, $`\varphi (𝐋)\varphi ^{}(𝐋^{})=𝒞_L^{\varphi \varphi }\delta (𝐋𝐋^{})/[2\pi L(L+1)]`$, which can be inferred from the temperature and polarization map 4-point functions hu02b . In Figure 1 we plot the deflection angle power spectrum, $`𝒞_l^{dd}l(l+1)𝒞_l^{\varphi \varphi }`$.

We calculate the 2-point functions using a publicly available code, CMBfast seljak96 , which was modified to include a scalar field dark energy component, to calculate $`𝒞_l^{dd}`$, and to include the effect of massive neutrinos on the recombination history (through the expansion rate). We use the Peacock and Dodds prescription to calculate the non-linear matter power spectrum peacock96 .

Effect of neutrinos. The lower panel in Figure 1 shows the differences in the power spectra between our fiducial model and the exact same model but with one of the three neutrino masses altered from zero to 0.1 eV. The error boxes are those for CMBpol (described below; see Table 1). The $`𝒞_l^{dd}`$ are noise-dominated at $`l>600`$ for CMBpol.

The signature of a 0.1 eV neutrino in the angular power spectra, in the absence of lensing, is at the 0.1% level. Such small masses are only detectable through their effect on lensing, which comes through their influence on the gravitational potential. Replacing a massless component with a massive one increases the energy density and therefore the expansion rate, suppressing growth. The net suppression of the power spectrum is scale dependent and the relevant length scale is the Jeans length for neutrinos bond83 ; ma96 ; hu98  which decreases with time as the neutrino thermal velocity decreases. This suppression of growth is ameliorated at scales larger than the Jeans length at matter–radiation equality, where the neutrinos can cluster. Neutrinos never cluster at scales smaller than the Jeans length today. The net result is no effect on large scales and a suppression of power on small scales, resulting in the shape of $`\delta 𝒞_l^{dd}/𝒞_l^{dd}`$ in Figure 1.

Error forecasting method.

The power spectra we include in our analysis are $`\stackrel{~}{𝒞}_l^{TT}`$, $`\stackrel{~}{𝒞}_l^{TE}`$, $`\stackrel{~}{𝒞}_l^{EE}`$ (unlensed), and $`𝒞_l^{dd}`$. We do not use the lensed power spectra to avoid the complication of the correlation in their errors between different $`\mathrm{}`$ values and with the error in $`𝒞_l^{dd}`$. Using the lensed spectra and neglecting these correlations can lead to overly optimistic forecasts hu02a . If we include the lensed spectra instead of the unlensed ones, the expected errors on $`w_x`$ and $`m_\nu `$ for CMBpol (see Table 1) shrink by about 40% and 30% respectively.

The distortions to the angular power spectra due to a $`0.1`$ eV neutrino and changes of order 10% in $`w_x`$ are very small. We have taken care to accurately forecast the constraints possible in this mass range. First, we make a Taylor expansion of the power spectra to first order in all the cosmological parameters. Then, given the the expected experimental errors on the power spectra, the expected parameter error covariance matrix is easily calculated.

The Taylor expansion works better and susceptibility to numerical error is reduced with a careful choice of the parameters used to span a given model space eisenstein99 ; efstathiou99 ; hu01a ; kosowsky02 . We take our set to be $`𝒫=\{\omega _m,\omega _b,\omega _\nu ,\theta _s,w_x,z_{\mathrm{ri}},k^3P_\mathrm{\Phi }^i(k_f),n_s,n_s^{},y_{\mathrm{He}}\}`$, with the assumption a flat universe. The first three of these are the densities today (in units of $`1.88\times 10^{29}\mathrm{g}/\mathrm{cm}^3`$) of cold dark matter plus baryons, baryons and massive neutrinos. Next two are the angular size subtended by the sound horizon on the last–scattering surface and the ratio of dark energy pressure to density. The Thompson scattering optical depth for CMB photons, $`\tau `$, is parameterized by the redshift of reionzation $`z_{\mathrm{ri}}`$. The primordial potential power spectrum is assumed to be $`k^3P_\mathrm{\Phi }^i(k)=k_f^3P_\mathrm{\Phi }^i(k_f)(k/k_f)^{n_s1+n_S^{}ln(k/k_f)}`$ with $`k_f=0.05\mathrm{Mpc}^1`$. The fraction of baryonic mass in Helium is $`y_{\mathrm{He}}`$. We Taylor expand about $`𝒫=\{0.146,0.021,0,0.6,1,6.3,6.4\times 10^{11},1,0,0.24\}`$.

We follow zaldarriaga97  to calculate the errors expected in $`\stackrel{~}{𝒞}_l^{TT}`$, $`\stackrel{~}{𝒞}_l^{TE}`$ and $`\stackrel{~}{𝒞}_l^{EE}`$ given Table 1. For errors on $`𝒞_l^{dd}`$ we follow hu02b . The errors on the unlensed spectra in the regime where lensing is important (deep in the damping tail) are certainly underestimated because reconstruction of the unlensed map from the lensed map will add to the errors. However, this is not worrisome since we limit all the unlensed spectra to $`l<2000`$, and a further restriction to $`l<1500`$ (where lensing is least important) only increases the error on $`m_\nu `$ by about 10% for CMBpol.

Experiments. We consider Planck planck , a high-resolution version of CMBpol 222CMBpol: http://spacescience.nasa.gov/missions/concepts.htm., and a polarized bolometer array on the South Pole Telescope 333SPT: http://astro.uchicago.edu/spt/ we will call SPTpol. Their specifications are given in Table 1. We assume that other frequency channels of Planck and CMBpol (not shown in the table) will clean out non-CMB sources of radiation perfectly. Detailed studies have shown foreground degradation of the results expected from Planck to be mild knox99 ; tegmark00b ; bouchet99 . At $`l>3000`$ emission from dusty galaxies will be a significant source of contamination. The effect is expected to be more severe for temperature maps. Hence we restrict temperature data to $`l<2000`$ and polarization data to $`l<2500`$.

Results. We emphasize the ability of the experiments to simultaneously determine $`P_\mathrm{\Phi }^i`$, $`w_x`$ and $`m_\nu `$ 444The degeneracy-breaking power of CMB lensing was first pointed out for $`\mathrm{\Omega }_\mathrm{\Lambda }/\mathrm{\Omega }_k`$ in metcalf98 .. These all affect the amplitude of $`P_\mathrm{\Phi }`$ at late times, the latter two due to their effect on the rate of growth of density perturbations. If we were only sensitive to the amplitude of $`𝒞_l^{dd}`$ then there would be an exact degeneracy between these three parameters. However, the $`l`$-dependence of the response of $`𝒞_l^{dd}`$ to these parameter variations breaks this would-be degeneracy, allowing for their simultaneous determination.

The effect of $`m_\nu `$ can easily be disentangled from that of $`w_x`$. We have already discussed the $`l`$-dependence of $`\mathrm{ln}𝒞_l^{dd}/m_\nu `$ shown in Fig. 1 as resulting from the scale- and time-dependence of $`\mathrm{ln}P_\mathrm{\Phi }/m_\nu `$. The $`l`$-dependence of $`\mathrm{ln}𝒞_l^{dd}/w_x`$ has the opposite sense. Although the suppression of $`P_\mathrm{\Phi }`$ for increasing $`w_x`$ is nearly $`k`$–independent, the effect is larger at late times —hence the radial projection gives a larger effect at low $`\mathrm{}`$. The effects of $`m_\nu `$ and $`w_x`$ are sufficiently distinct to allow for their simultaneous determination. We point out that the effect of $`w_x`$ is more pronounced for larger values due to two reasons. One, dark energy starts to dominate earlier (which implies larger uniform suppression) and two, perturbations in dark energy on large scales are enhanced for large $`w_x`$.

The difference in the response of $`𝒞_l^{dd}`$ to $`m_\nu `$ and $`w_x`$ allows for, e.g., Planck to detect the acceleration of the Universe ($`w_x<1/(3\mathrm{\Omega }_x)`$) at the 2$`\sigma `$ level. Such a confirmation would be valuable given the deep theoretical implications of acceleration kaloper02 . Hu hu02a  has previously noted this result obtained with the assumption $`m_\nu =0`$.

As is well known, the $`P_\mathrm{\Phi }^i`$ can be determined independently of the lensing signal, through use of a signal at large angular scales. One combines $`𝒞_l^{EE}`$ and $`𝒞_l^{TE}`$ at $`l\begin{array}{c}<\hfill \\ \hfill \end{array}20`$ where they are proportional to $`P_\mathrm{\Phi }^i\tau ^2`$ and $`P_\mathrm{\Phi }^i\tau `$ respectively zaldarriaga97 ; kaplinghat03a  with the TT, EE and ET spectra at $`20\begin{array}{c}<\hfill \\ \hfill \end{array}l\begin{array}{c}<\hfill \\ \hfill \end{array}2000`$ where they are proportional to $`P_\mathrm{\Phi }^ie^{2\tau }`$.

If we assume a single-step transition for the ionization history Planck can achieve $`\sigma (\tau )=0.005`$ eisenstein99 . However, foreground contamination tegmark00b , and modeling uncertainty in the ionization history holder03  can increase this uncertainty. For these reasons we conservatively ignore polarization data at $`l<30`$ and instead set a prior, by hand, of $`\sigma (\tau )=0.009`$; including the $`l<30`$ polarization data would (perhaps artificially) achieve a smaller $`\sigma (\tau )`$. In the end, $`\tau `$ is determined (only slightly) better than this prior because there is some constraint on $`P_\mathrm{\Phi }^i`$ from the lensing signal. Note that since $`P_\mathrm{\Phi }^ie^{2\tau }`$ is so well-determined, we always expect $`\sigma (\mathrm{ln}P_\mathrm{\Phi })=2\sigma (\tau )`$, as we find..

An extended period of reionization, as suggested by the combination of WMAP and quasar observations kogut03 , may have large spatial fluctuations in the ionization fraction. Such “patchy” reionization would lead to a large diffuse kinetic SZ contribution to $`𝒞_l^{TT}`$ at high $`l`$ knox98 ; santos03 , possibly larger than the lensing contribution. Fortunately the analogous effect in the polarization is much smaller. For a conservative upper bound on how patchy reionization could degrade $`\sigma (m_\nu )`$, we restrict the temperature data to $`l<1000`$ and find $`\sigma (m_\nu )=0.045`$ eV for CMBpol and 0.34 eV for Planck.

The primary motivation for CMBpol is the detection of the B mode due to gravity waves produced in inflation. The amplitude of this signal would directly tell us the energy density during inflation. Following the calculation in knox02 ; kesden02  we find a $`3\sigma `$ detection is possible for CMBpol if the energy density during inflation is greater than $`\rho _{\mathrm{min}}=(2\times 10^{15}\mathrm{GeV})^4`$; $`\rho _{\mathrm{min}}^{1/4}`$ is an order of magnitude smaller than the GUT scale. We note that $`\rho _{\mathrm{min}}1/\tau `$, approximately, for $`0.05<\tau <0.2`$ and we have assumed $`\tau =0.1`$. This scaling with $`\tau `$ suggests that the reionization feature in the B mode at the largest angular scales is important and therefore a full-sky experiment is necessary to achieve this sensitivity level.

The scalar spectrum determined from high–resolution CMB observations (the constraining power comes from primary CMB) can also be a useful probe of inflation, as studied recently by gold03 . If $`n_S1=0.07`$, the central value in fits to WMAP and other observations spergel03 , then inflationary models generically predict $`n_S^{}(n_S1)^2=0.005`$ which will be detectable at the $`3\sigma `$ level by CMBpol.

Determining $`\omega _b`$ and $`y_{He}`$ to high precision will facilitate precision consistency tests with Big Bang Nucleosynthesis (BBN) predictions. It will also be useful in constraining non-standard BBN. For example, determining $`\omega _b`$ and $`y_{He}`$ to high precision allows strong constraints to be put on the number of relativistic species N (or equivalently the expansion rate) during BBN. If $`\sigma (y_{He})`$ is small, then $`\sigma (N)=\sigma (y_{He})/0.013`$, which for CMBpol works out to $`\sigma (N)=0.4`$. Constraints on $`N`$ have important repercussions for neutrino mixing in the early universe, and hence on neutrino mass models abazajian02 .

Conclusions. Gravitational lensing of the CMB is a promising probe of the growth of structure and the fundamental physics that affects it. High sensitivity, high resolution maps will allow us to measure the lensing signature well enough to simultanteously constrain $`m_\nu `$, $`w_x`$ and $`P_\mathrm{\Phi }`$. A future all-sky polarized CMB mission aimed at detecting gravitational waves is likely to succeed in determining neutrino mass as well.

###### Acknowledgements.

We thank J. Bock, S. Church, W. Hu, M. Kamionkowski, A. Lange, S. Meyer and M. White for useful conversations.


text:
# Numerical analysis of a time-headway bus route model

## I Introduction

While there has been much interest in the study of automobile traffic Traffic , there have been few corresponding studies of buses Buses ; Nagatani-old ; Nagatani ; Nagatani-new . The dynamics of a bus route, while having some similarities with that of automobile traffic, differs due to the added interaction of buses with passengers at designated bus stops. A good reason for studying the dynamics of bus routes is that they are so often unstable. Buses are initially spaced at regular intervals. However, if one bus is delayed for some reason, it will find a larger number of passengers waiting for it at subsequent stops, delaying it further. Meanwhile, the bus following finds fewer passengers waiting for it, allowing it to go faster until eventually it meets up with the delayed bus. Clusters of three, four, or more buses have been known to form in this manner, resulting in slower service.

In references Nagatani-old  and Nagatani , Nagatani presents a time-headway model for buses. Using linear stability analysis, he is able to determine the range of parameters over which the homogeneous solution (i.e., with buses spaced evenly apart) is unstable. In this paper, we make a more thorough investigation of Nagatani’s model. We demonstrate the existence of three types of phase diagrams, in which the behavior of the bus system is divided into four separate categories. We conclude with a discussion of how this model may be applied to real-world situations, and the limitations imposed by practical considerations.

## II Model

We consider the following model, introduced in Nagatani , of buses on a bus route (Fig. 1). Bus stops are labelled by $`s=1,2,\mathrm{}`$ where stops $`s`$ and $`s+1`$ are a distance $`L`$ apart. There are $`J`$ buses, $`j=1,\mathrm{},J`$, which travel from stop to stop, with bus $`j=1`$ in the lead and bus $`j=J`$ in the rear. Every bus visits every stop, and buses do not pass one another. The time at which bus $`j`$ arrives at stop $`s`$ is $`t_{j,s}`$, which is given by the recursive relation

$$t_{j,s}t_{j,s1}=\lambda \gamma \delta t_{j,s1}+\frac{L}{V_{j,s1}},$$
(1)
where

$$\delta t_{j,s}t_{j,s}t_{j1,s}$$
(2)
is the *time-headway*, the time gap in front of bus $`j`$ at stop $`s`$. The first term on the right-hand side of Eq. 1 is the time it takes for passengers to board the bus at stop $`s1`$. The parameter $`\lambda `$ is the rate at which passengers arrive at a bus stop; $`\lambda \delta t_{j,s1}`$ is the number of passengers that have arrived at stop $`s1`$ since the previous bus left. The parameter $`\gamma `$ is the time it takes each passenger to board the bus, so $`\lambda \gamma \delta t_{j,s1}`$ is the amount of time needed to board all of the passengers. For convenience, we introduce the dimensionless parameter $`\mu \lambda \gamma `$, which we call the *passenger rate*. For simplicity we ignore the passengers getting off of the bus; we will assume that it takes less time for the passengers to get off than it does to get on and pay the fare.

The second term in Eq. 1 is the time it takes for bus $`j`$ to travel from stop $`s1`$ to stop $`s`$, where $`V_{j,s1}`$ is the average velocity of the bus between stops. If this velocity is constant, then the tendency for buses to bunch together, as described in the introduction, has no counterweight, and a steady flow of buses will always be unstable (unless there are no passengers). It is reasonable to assume, however, that a bus driver will try to prevent bunching by slowing down when the gap between his bus and the next is too small. One can model this by writing the average speed $`V_{j,s}`$ as a function $`\stackrel{~}{V}(\delta t_{j,s})`$ of the gap between his bus and the bus in front of him:

$$\stackrel{~}{V}(\delta t)=v_{\mathrm{min}}+(v_{\mathrm{max}}v_{\mathrm{min}})\frac{\mathrm{tanh}\omega (\delta tt_c)+\mathrm{tanh}\omega t_c}{1+\mathrm{tanh}\omega t_c}$$
(3)
The hyperbolic tangent factor acts as a spread-out step function, centered at $`t_c`$ with a width proportional to $`1/\omega `$. The parameter $`v_{\mathrm{max}}`$ is the speed a free bus (i.e.one that is alone on the route) would travel. On the other hand, $`v_{\mathrm{min}}`$ is the speed a bus travels if it has completely caught up with the bus in front of it. For example, if $`v_{\mathrm{min}}=0`$ then a bus which has caught up with the bus in front of it will stop and wait until its predecessor has cleared the next stop before proceeding.

In what follows, it is convenient to work with the time headways $`\delta t_{j,s}`$, rather than the arrival times $`t_{j,s}`$. It is also convenient to rewrite our expressions in terms of dimensionless quantities. In doing so, we find that there are four significant parameters, not including initial conditions. The first such parameter is the passenger rate $`\mu `$. The other three are $`\alpha =L\omega /v_{\mathrm{max}}`$, $`\beta =v_{\mathrm{min}}/v_{\mathrm{max}}`$, and $`ϵ=1\mathrm{tanh}\omega t_c`$ (which will typically be small). We will also consider the dimensionless variable $`\mathrm{\Delta }t_{j,s}=\omega \delta t_{j,s}`$ and the dimensionless velocity function $`V(\mathrm{\Delta }t)=\stackrel{~}{V}(\mathrm{\Delta }t)/v_{\mathrm{max}}`$. Our evolution equation Eq. 1 now reads

$$\mathrm{\Delta }t_{j,s}\mathrm{\Delta }t_{j,s1}=\alpha \left[\frac{1}{V(\mathrm{\Delta }t_{j,s1})}\frac{1}{V(\mathrm{\Delta }t_{j1,s1})}\right]+\mu [\mathrm{\Delta }t_{j,s1}\mathrm{\Delta }t_{j1,s1}],$$
(4)
where

$$V(\mathrm{\Delta }t)=\beta +\frac{(1\beta )ϵ\mathrm{tanh}\mathrm{\Delta }t}{1(1ϵ)\mathrm{tanh}\mathrm{\Delta }t}=\frac{\beta (1\mathrm{tanh}\mathrm{\Delta }t)+ϵ\mathrm{tanh}\mathrm{\Delta }t}{(1\mathrm{tanh}\mathrm{\Delta }t)+ϵ\mathrm{tanh}\mathrm{\Delta }t}.$$
(5)
## III Stability Analysis

We are interested in the stability of a homogeneous flow of buses, with $`\mathrm{\Delta }t_{j,s}=\mathrm{\Delta }t_{j,0}=\mathrm{\Delta }t_0`$. One can easily verify that this is a solution to Eq. 4. Starting with a small perturbation to the initial homogeneous solution: $`\mathrm{\Delta }t_{j,s}=\mathrm{\Delta }t_0+y_{j,s}`$, where $`y_{j,s}`$ is small. To first order, Eq. 4 becomes

$$y_{j,s}y_{j,s1}=[y_{j+1,s1}y_{j,s1}][F(\mathrm{\Delta }t_0)\mu ],$$
(6)
where we have introduced the convenient abbreviation

$$F(\mathrm{\Delta }t_0)\alpha \frac{V^{}(\mathrm{\Delta }t_0)}{V(\mathrm{\Delta }t_0)^2}=\frac{\alpha (1\beta )ϵ(1\mathrm{tanh}^2\mathrm{\Delta }t_0)}{[\beta (1\mathrm{tanh}\mathrm{\Delta }t_0)+ϵ\mathrm{tanh}\mathrm{\Delta }t_0]^2}.$$
(7)
It can be shown Nagatani  that the perturbation is stable if

$$F(\mathrm{\Delta }t_0)1<\mu <F(\mathrm{\Delta }t_0).$$
(8)
From Eq. 8 we can construct a phase diagram (Fig. 2) for the stability of an initially homogeneous bus route, based on the initial spacing $`\mathrm{\Delta }t_0`$ and the passenger rate $`\mu `$. The stable region in phase space is bounded by the curves $`\mu =F(\mathrm{\Delta }t_0)`$ and $`\mu =F(\mathrm{\Delta }t_0)1`$. Because of the added constraint that $`\mu 0`$, there are different phase diagrams depending on whether $`F(\mathrm{\Delta }t_0)1`$ is ever positive (Fig. 2a) or not (Fig. 2b). The curve $`F(\mathrm{\Delta }t)`$ has a maximum value of

$$\alpha \frac{(1\beta )}{2\beta ϵ}\text{ at }x=1\frac{ϵ}{\beta }.,$$
(9)
so the phase diagram resembles Fig. 2a whenever

$$\alpha >\frac{2\beta ϵ}{1\beta }\frac{2\beta }{1\beta }.$$
(10)
A third phase diagram, Fig. 2c, arises when $`v_{\mathrm{min}}=0`$, as it is in figure 3 of Ref. Nagatani  (although apparently not in figure 8 of the same reference, which may account for the discrepancy between those two phase diagrams.)

## IV Simulation

To study the ways in which the system becomes unstable, we evaluate Eq. 4 iteratively in $`s`$. Our initial condition is

$$\mathrm{\Delta }t_{j,0}=\mathrm{\Delta }t_0+0.1r_j,$$
(11)
where $`r_j`$ are random numbers chosen between $`1`$ and $`1`$. For each combination of initial headway $`\mathrm{\Delta }t_0`$ and passenger rate $`\mu `$, we run the simulation until either a) we reach stop $`s=5000`$, or b) one or more of the bus headways exceeds $`\mathrm{\Delta }t=1000`$ (in which case the system has become unphysical).

In this paper we consider two different boundary conditions. The first is periodic in the bus number $`j`$; so for example $`\mathrm{\Delta }t_{1,s}=t_{1,s}t_{J,s}`$. This is convenient numerically, and it creates translational symmetry, but it is hard to construct a physical model which begins with this characteristic. We also consider a fixed boundary condition, where $`\mathrm{\Delta }t_{1,s}=\mathrm{\Delta }t_0`$. Since the velocity of a bus depends entirely on $`\mathrm{\Delta }t`$, this corresponds to a scenario where the initial bus ($`j=1`$) moves at a constant speed $`V(\mathrm{\Delta }t_0)`$.111 A more realistic boundary condition may involve ramping up the passenger rate $`\mu `$ over time, as occurs during the course of a normal day.

The structure of the model requires that buses not pass one another; however, there is nothing in Eq. 4 to prevent the headways $`\mathrm{\Delta }t`$ from becoming negative. To fix this, we add to our simulation the rule that any $`\mathrm{\Delta }t_{j,s}<0`$ is replaced by $`\mathrm{\Delta }t_{j,s}=0`$. This corresponds to a situation where drivers are forbidden (or unable due to road conditions) to pass one another.222An alternative solution which we do not consider here is to allow buses to pass one another. This could conceivably be done by replacing all $`\mathrm{\Delta }t_j<0`$ with $`\left|\mathrm{\Delta }t_j\right|`$, effectively swapping the labels of buses that pass one another. We have not considered whether this would work in practice, however.

Figure 3 shows the results of our simulation runs for a typical set of parameters ($`\alpha =1`$, $`\beta =1/4`$, $`ϵ=1\mathrm{tanh}2=0.036`$), using both types of boundary conditions. In both cases, the phase space is divided into four regions, corresponding to four types of runs.

### IV.1 Stable Runs

Most of the runs within the stable region, as defined by Eq. 8, remain homogeneous. In the periodic case, the initial fluctuations in $`\mathrm{\Delta }t_j`$ settle into a small precessing sinusoidal perturbation which decays exponentially with time (Fig. 4). In the fixed case, the system quickly locks onto the constant solution $`\mathrm{\Delta }t=\mathrm{\Delta }t_0`$ with no fluctuations.

### IV.2 Explosive Runs

Most of the runs lying above the stable region quickly develop an unphysical instability. This takes the form seen in Fig. 5, independent of boundary condition: those headways lying above the mean increase exponentially, while those lying below decrease steadily until they reach zero. An observer stationed at a stop far down the line will see clusters of buses arriving after long waits; far enough down the line, these waits become astronomical, which is absurd. Clearly this model is insufficient to deal with these runs at long times.

### IV.3 Slowed Runs

To the left of the stable region are runs which develop an alternative stable solution, as seen in Fig. 6. In the case of the fixed boundary condition, these runs have two things in common. The first, indicated by the vanishing of one or more headways, is the appearance of clusters: two or more buses which travel along as a single unit. The second is that the units, whether single buses or clusters, are homogeneously spaced, but with a spacing that is larger than the initial spacing $`\mathrm{\Delta }t_0`$. It should be pointed out that the solutions shown here are stationary; the clusters and spacings, after they form, do not change.

In analytic terms, these states are of the form $`\mathrm{\Delta }t_{j,s}=\mathrm{\Delta }t_j=\tau r_j`$, where $`r_j`$ is either $`0`$ or $`1`$, and $`\tau >\mathrm{\Delta }t_0`$ is a constant. It is straightforward to show that this is a solution to Eq. 4:

$$0=\alpha \left[\frac{1}{V(\tau r_j)}\frac{1}{V(\tau r_{j1})}\right]+\mu \tau [r_jr_{j1}]$$
(12)
When $`r_j=r_{j1}`$, this equation is satisfied trivially. Otherwise, the equation takes the form

$$\mu =\frac{\alpha }{\tau }\left(\frac{1}{\beta }\frac{1}{V(\tau )}\right).$$
(13)
which we can solve numerically for $`\tau `$ (Fig. 7). For a given passenger rate $`\mu `$, these spacings $`\tau `$ correspond precisely with those seen in simulation. Furthermore, for high enough passenger rates— $`\mu >1.199`$ for this set of parameters— Eq. 13 has no real solutions, which explains the cut-off in Fig. 3 between the slowed and explosive regimes.

In the case of the periodic boundary condition there are cases where the clusters eventually break up, leaving a system of buses which are equally spaced, but with the larger spacing predicted by Eq. 13. These runs are marked by downward-pointing triangles ($``$) in Fig. 3.

### IV.4 Oscillatory Runs

In the case of a run lying below the stability region in phase space, the first term in Eq. 4, which is meant to resist the tendency for buses to cluster, becomes too large. This leads to overreaction, so that two buses which arrive too close together at one stop are too far apart at the next. The resulting behavior may be compared to a system of underdamped oscillators. Figure 8 shows the resulting behavior. For the periodic boundary condition, these oscillations decay as a power law(Fig. 9), but at so slow a rate as to be practically permanent.

With the fixed boundary condition, the earliest buses (i.e., those with the lowest $`j`$) shed the oscillating behavior after only a few stops, resuming a homogeneous configuration; with more iterations, more buses join the homogeneous regime. In some cases, such as in Fig. 10, the system reaches a steady state with the oscillations still dominating the later buses. In other cases, however, the system becomes completely homogeneous, as the bifurcation point seen in Fig. 10 slips off the right side of the graph. This later effect may be due to the finite number of buses.

## V Discussion

In our simulations we have considered a large range of values for $`\mu `$ and $`\mathrm{\Delta }t_0`$. However, these parameters should be limited by a couple of practical concerns.

The passenger rate $`\mu `$ is defined as the product of the number of passengers that arrive at a stop per unit time, and the time it takes a single person to board the bus. Put another way, it is the ratio of the number of people that arrive at a stop to the number of people who can board the bus in the same amount of time. If this number is greater than $`1`$, then passengers arrive at a stop faster than the bus can take them on, and the bus should never be able to leave the stop. Since $`0<\mu <1`$, only one of the two inequalities in Eq. 8 is meaningful for a given value of the parameters (including $`\mathrm{\Delta }t_0`$), depending on whether $`F(\mathrm{\Delta }t_0)`$ is larger or smaller than 1. This suggests that if one wanted to maximize the area of the stability region in phase space, one would do well to make sure that the lower stability curve $`F(\mathrm{\Delta }t_0)1`$ just grazes zero, or that $`\alpha (1\beta )=2\beta ϵ`$ according to Eq. 9.

Another practical consideration puts a limit on the value of $`\mathrm{\Delta }t_0`$. Typically, buses are spaced far enough apart so that the first bus will reach the first stop before the second bus is allowed to leave, particularly if the stops are spaced fairly close together. This is described by the inequality

$$\mathrm{\Delta }t_0>\frac{L}{\stackrel{~}{V}(\mathrm{\Delta }t_0)}=\alpha \frac{(1\mathrm{tanh}\mathrm{\Delta }t_0)+ϵ\mathrm{tanh}\mathrm{\Delta }t_0}{\beta (1\mathrm{tanh}\mathrm{\Delta }t_0)+ϵ\mathrm{tanh}\mathrm{\Delta }t_0}>\alpha .$$
(14)
For the parameters we have been studying,

$$\mathrm{\Delta }t_0>\frac{10.964\mathrm{tanh}\mathrm{\Delta }t_0}{0.250.214\mathrm{tanh}\mathrm{\Delta }t_0}\mathrm{\Delta }t_0>1.82.$$
(15)
This cuts out much of the interesting part of Fig 3, including the slowed runs and almost all of the underdamped solutions. In our discrete model, the basic iteration step is the bus stop; bus drivers are allowed to change their speed at the bus stops and nowhere else. If the buses are several stops apart, then they have enough time to react to one another. Otherwise, unusual situations such as the slowing case or the underdamped case may arise.

Finally, we consider the relationship between $`\delta t_0`$ and $`t_c`$, which is how close a bus will come to the bus in front of it without slowing down. Consider the ratio

$$r=\frac{1\mathrm{tanh}\mathrm{\Delta }t_0}{ϵ}=\frac{1\mathrm{tanh}\omega \delta t_0}{1\mathrm{tanh}\omega t_c}.$$
(16)
We can rewrite $`F(\mathrm{\Delta }t_0)`$ in terms of this ratio:

$$F(\mathrm{\Delta }t_0)=\alpha (1\beta )\frac{r(2ϵr)}{[\beta r+1ϵr]^2}.$$
(17)
If $`ϵr=1\mathrm{tanh}\mathrm{\Delta }t_01`$ (which it will be if Eq. 15 is valid, since $`1\mathrm{tanh}1.82=0.05`$), then

$$F(\mathrm{\Delta }t_0)\alpha (1\beta )\frac{2r}{(1+\beta r)^2}.$$
(18)
Now let us consider what values our parameters might take in real life. A typical urban bus route might have $`L=0.5\mathrm{km}`$, $`v_{\mathrm{max}}=50\mathrm{km}/\mathrm{hr}=\frac{5}{6}\mathrm{km}/\mathrm{min}`$, $`v_{\mathrm{min}}=15\mathrm{km}/\mathrm{hr}=\frac{1}{4}\mathrm{km}/\mathrm{min}`$, and $`\omega =1/\mathrm{min}`$; thus $`\alpha =0.6`$ and $`\beta =0.3`$. A bus which runs every 10 minutes might take on two passengers at every stop, so $`\lambda =0.2`$ people per minute. If it takes $`\gamma =3\mathrm{s}`$ for a person to board a bus, then $`\mu =\lambda \gamma =0.01`$.

Consider a scenario where bus drivers only react to what they see; that is, a driver will only slow down if she can see the next bus in front of her. It takes a free bus $`L/v_{\mathrm{max}}=0.6`$ minutes to travel from one stop to the next, so a reasonable value for the amount of warning a bus driver has is on the order of $`t_c=1\mathrm{min}`$. Typical bus routes have buses which are spaced much farther apart, perhaps every $`\mathrm{\Delta }t_0=10\mathrm{min}`$ or more. In this scenario, $`r=(1\mathrm{tanh}10)/(1\mathrm{tanh}1)=10^8`$, so $`F(\mathrm{\Delta }t_0)10^8`$. Since $`\mu =10^2`$, the stability condition Eq. 8 is very clearly violated. For $`F(\mathrm{\Delta }t_0)`$ to reach a high enough value to create stability, we need in general for the ratio $`r`$ to be closer to $`1`$. $`F(r)`$ takes its maximum value when $`r=1/(2\beta )`$, in which case $`F(r)=1/(8\beta )=0.4`$, which is easily larger than $`\mu `$ in this example. Notice that, for values of $`t_c`$ and $`\delta t_0`$ greater than $`1\mathrm{min}`$, $`re^{2(t_c\delta t_0)}`$, so for each minute’s difference between $`t_c`$ and $`\delta t_0`$, $`r`$ is increased or decreased by a factor of ten. It would seem that to maintain a stable homogeneous bus route, drivers must be reacting to the bus in front of them even from the very beginning, and can only ignore the leading bus if they have gotten far enough behind (in which case, of course, the proper solution is to go as quickly as possible).

This is quite a stringent requirement for stability, and explains why it is so common to see clusters of buses in large cities. It does not seem likely, however, that this would be the case for less frequent bus service, such as when buses run once per hour. A driver on such a route does not typically keep track of what the previous bus was doing an hour ago, and yet one does not see clustering behavior on these low-frequency routes. The reason for this is that the instabilities predicted by Eq. 8 may take a long time to become noticeable, and normal bus routes tend to have a limited number of stops. Figure 11 shows the number of stops that a bus route has to cover before seeing a noticeable deviation in the initial homogeneous state. If the passenger rate is $`\mu =0.01`$ as suggested in the urban case above, then the route would have to have 130 stops to show a 1 minute deviation from the homogeneous state, and 225 stops to show a 5 minute deviation. If $`L=0.5\mathrm{km}`$, these correspond to $`65\mathrm{km}`$ and $`112\mathrm{km}`$, longer than your average bus route. The situation is even better when you consider that a suburban or rural bus route might have, not 1 person for every 5 minutes, but maybe 1 person every 25 minutes, so that $`\mu =0.002`$, and we can start having bus routes with $`600`$ or $`1000`$ stops before the instability becomes noticeable. This is not to say that smaller routes remain perfectly on time, of course; just that the delays are unlikely to be due to the need to pick up extra passengers. Since buses will typically complete the route only to turn around and do it again, one might consider an entire day’s run to be a single route, in which case instabilities may creep in late in the day. However, the introduction of a bus terminal, where buses wait until a specific time to leave for their next trip through the route, would have to be accounted for in this case.

In this paper, we have considered the bus route model proposed in Ref. Nagatani . We have added a simple way to deal with negative time-headways (by replacing all negative $`\mathrm{\Delta }t_j`$’s with zeroes), and by doing so have been able to study the unstable modes of a homogeneous system of buses. We show that there are in fact three different phase diagrams (Fig 2) for this system, depending on our choice of parameters, and that in addition to the stable homogeneous state, there are three unstable modes which the system can fall into: the explosive mode, the slowed mode, and the underdamped mode. We then considered the application of this model to real-life bus routes. We have shown that the passenger rate $`\mu `$ and the initial spacing between buses $`\mathrm{\Delta }t_0`$ are greatly restricted by practical considerations, and that to guarantee stability one needs to have bus drivers who are constantly tracking the bus in front of them, even when that bus is at its normal distance. Fortunately, this is only necessary for bus routes with very many stops; with fewer stops to make, a bus may be able to complete the route before any instabilities can become noticeable.

The author thanks Dr. Raghuveer Parthasarathy for indirectly and unknowingly encouraging this work; and Professor Gene Mazenko for the computing power needed to run the simulations. This work was supported by the Materials Research Science and Engineering Center through Grant No. NSF DMR 9808595.


text:
# An Infrared Renormalization Group Limit Cycle in QCD

(March 17, 2003)

## Abstract

Small increases in the up and down quark masses of QCD would tune the theory to the critical renormalization group trajectory for an infrared limit cycle in the three-nucleon system. At critical values of the quark masses, the deuteron binding energy goes to zero and the triton has infinitely many excited states with an accumulation point at the 3-nucleon threshold. The ratio of the binding energies of successive states approaches a universal constant that is close to 515. The proximity of physical QCD to the critical trajectory for this limit cycle explains the success of an effective field theory of nucleons with contact interactions only in describing the low-energy 3-nucleon system.

Renormalization group, limit cycle, quantum chromodynamics

The development of the renormalization group (RG) has had a profound effect on many branches of physics. Its successes range from explaining the universality of critical phenomena in condensed matter physics to the nonperturbative formulation of quantum field theories that describe elementary particles Wilson:dy . Most simple applications of the renormalization group involve fixed points that are approached either in the infrared limit or in the ultraviolet limit. However, the renormalization group can be reduced to a set of differential equations that define a flow in the space of coupling constants. A fixed point where the flow vanishes is only the simplest topological feature that can be exhibited by such a flow. Another topological feature is a limit cycle, which is a 1-dimensional orbit that is closed under the RG flow. The possibility of RG flow to a limit cycle was proposed by Wilson in 1971 Wilson:1970ag . Glazek and Wilson have recently constructed a simple discrete Hamiltonian system that exhibits RG flow to a limit cycle Glazek:2002hq . However, few physical applications of RG limit cycles have emerged.

The purpose of this letter is to point out that Quantum Chromodynamics (QCD) is close to the critical trajectory for an infrared RG limit cycle in the 3-nucleon sector. It can be tuned to the critical trajectory by small changes in the up and down quark masses. The proximity of the physical quark masses to these critical values explains the success of a program initiated by Efimov to describe the 3-nucleon problem in terms of zero-range forces between nucleons Efi81 . An effective-field-theory formulation of this program by Bedaque, Hammer, and van Kolck exhibits an ultraviolet RG limit cycle Bedaque:1999ve . The proximity of physical QCD to the critical trajectory for an infrared limit cycle shows that the ultraviolet limit cycle of Ref. Bedaque:1999ve  is not just an artifact of their model.

In the late 1960’s, Wilson used the renormalization group (RG) to explain universality in critical phenomena Wilson:dy . Transformations of a system that integrate out short-distance degrees of freedom while leaving the long-distances physics invariant define a RG flow on the multidimensional space of coupling constants $`𝐠`$ for all possible terms in the Hamiltonian:

$`\mathrm{\Lambda }{\displaystyle \frac{d}{d\mathrm{\Lambda }}}𝐠=𝜷(𝐠),`$ (1)
where $`\mathrm{\Lambda }`$ is an ultraviolet momentum cutoff. Standard critical phenomena are associated with infrared fixed points $`𝐠_{}`$ of the RG flow, which satisfy $`𝜷(𝐠_{})=0`$. The tuning of macroscopic variables to reach a critical point corresponds to the tuning of the coupling constants $`𝐠`$ to a critical trajectory that flows to the fixed point $`𝐠_{}`$ in the infrared limit $`\mathrm{\Lambda }0`$. One of the signatures of an RG fixed point is scale invariance: symmetry with respect to the coordinate transformation $`𝐫\lambda 𝐫`$ for any positive number $`\lambda `$. This symmetry implies that dimensionless variables scale as powers of the momentum scale, perhaps with anomalous dimensions.

In a classic 1971 paper, Wilson suggested that the renormalization group might also be relevant to the strong interactions of elementary particle physics Wilson:1970ag . At that time, the fundamental theory for the strong interactions had not yet been discovered. It was believed to involve quarks, and hints that the strong interactions might have scaling behavior at high energies had been observed in experiments on deeply-inelastic lepton-nucleon scattering. Wilson suggested that simple high-energy behavior can be explained by simple RG flow of the relevant coupling constants in the ultraviolet limit $`\mathrm{\Lambda }\mathrm{}`$. The simplest possibility is RG flow to an ultraviolet fixed point. Another simple possibility is RG flow to an ultraviolet limit cycle. A limit cycle is a 1-parameter family of coupling constants $`𝐠_{}(\theta )`$ that is closed under the RG flow and can be parametrized by an angle $`0<\theta <2\pi `$. The RG flow carries the system around a complete orbit of the limit cycle every time the ultraviolet cutoff $`\mathrm{\Lambda }`$ increases by some factor $`\lambda _0`$. One of the signatures of an RG limit cycle is discrete scale invariance: symmetry with respect to the coordinate transformation $`𝐫\lambda _0^n𝐫`$ only for integer values of $`n`$. This symmetry implies that dimensionless variables with vanishing anomalous dimensions are periodic functions of the logarithm of the momentum scale with period $`\mathrm{ln}(\lambda _0)`$. The fundamental field theory for the strong interactions, QCD, was eventually discovered. QCD has a single coupling constant $`\alpha _s(\mathrm{\Lambda })`$ with an asymptotically-free ultraviolet fixed point: $`\alpha _s(\mathrm{\Lambda })0`$ as $`\mathrm{\Lambda }\mathrm{}`$ Gross-Politzer . The renormalization of QCD does not involve an ultraviolet limit cycle.

We now turn to an independent theoretical development whose connection with RG limit cycles was first pointed out in Ref. Albeverio:zi . In 1970, Efimov discovered a remarkable effect that can occur in the 3-body sector for nonrelativistic particles with a resonant short-range 2-body interaction Efi71 . If the particles are identical bosons and if the resonance is in the S-wave channel, the strength of the interaction is governed by the S-wave scattering length $`a`$. Efimov showed that if $`|a|`$ is much larger than the range $`r_0`$ of the interaction, there are shallow 3-body bound states whose number increases logarithmically with $`|a|/r_0`$. In the resonant limit $`a\pm \mathrm{}`$, there are infinitely many shallow 3-body bound states with an accumulation point at the 3-body scattering threshold. The ratio of the binding energies of successive states rapidly approaches the universal constant $`\lambda _0^2515`$, where $`\lambda _0=e^{\pi /s_0}22.7`$ and $`s_01.00624`$ is a transcendental number. In subsequent work, Efimov studied both the bound state spectrum and low-energy 3-body scattering observables for identical bosons with large scattering length Efi79 . He showed that the observables for different values of $`a`$ are related by a discrete scaling transformation in which $`a\lambda _0^na`$, where $`n`$ is an integer, and lengths and energies are scaled by the appropriate powers of $`\lambda _0^n`$.

The Efimov effect can also occur for fermions with equal masses and large scattering lengths if they have at least 3 distinct spin states. Nucleons are examples of fermions with large scattering lengths. The spin-singlet and spin-triplet $`np`$ scattering lengths are $`a_s=23.8`$ fm and $`a_t=5.4`$ fm. They are both significantly larger than the $`NN`$ effective range, which is $`r_0=1.8`$ fm in the spin-triplet channel. Efimov used this observation as the basis for a qualitative approach to the 3-nucleon problem Efi81 . His starting point was to take nucleons as point particles with zero-range potentials whose strengths are adjusted to reproduce the scattering lengths $`a_s`$ and $`a_t`$. The effective range and higher order terms in the low-energy expansions of the phase shifts were treated as perturbations. This approach works well in the 2-nucleon system, giving accurate predictions for the deuteron binding energy. This is no surprise; it simply reflects the well-known success of the effective range expansion in the 2-nucleon system Bet49 . Remarkably, Efimov’s program also works well in the 3-nucleon system at momenta small compared to $`m_\pi `$. In the triton channel, the Efimov effect makes it necessary to impose a boundary condition on the wavefunction at short distances. The boundary condition can be fixed by using either the neutron-deuteron scattering length or the triton binding energy as input. But if one of them is used as input, the other is predicted with surprising accuracy. The accuracy can be improved by taking into account the effective range as a first-order perturbation Efi91 .

The Efimov effect was revisited by Bedaque, Hammer, and van Kolck within the framework of effective field theory (EFT) Bedaque:1998kgkm . The problem of bosons with mass $`m`$ and large scattering length $`a`$ can be described by a nonrelativistic field theory with a complex-valued field $`\psi `$ and Hamiltonian density

$`=(\mathrm{}^2/2m)\psi ^{}\psi +g_2(\mathrm{\Lambda })(\psi ^{}\psi )^2.`$ (2)
For convenience, we set $`\mathrm{}=1`$ in the following. In the 2-body sector, the nonperturbative solution of the field theory can be obtained analytically. Renormalization can be implemented by adjusting the 2-body coupling constant $`g_2(\mathrm{\Lambda })`$ as a function of the ultraviolet cutoff $`\mathrm{\Lambda }`$ so that the scattering length is $`a`$. Other 2-body observables are then independent of $`\mathrm{\Lambda }`$ and have the appropriate values for bosons with zero effective range.

In the 3-body sector, the nonperturbative solution of the field theory can be obtained by solving integral equations numerically. These integral equations have unique solutions only in the presence of an ultraviolet cutoff $`\mathrm{\Lambda }`$. The resulting predictions for 3-body observables, although finite, depend on the ultraviolet cutoff and are periodic functions of $`\mathrm{ln}(\mathrm{\Lambda })`$ with period $`\pi /s_0`$. Bedaque, Hammer, and van Kolck showed that the quantum field theory could be fully renormalized to remove the residual dependence on $`\mathrm{\Lambda }`$ in the 3-body sector by adding a 3-body interaction term $`g_3(\mathrm{\Lambda })(\psi ^{}\psi )^3`$ to the Hamiltonian density in (2) Bedaque:1998kgkm . The dependence of 3-body observables on the cutoff decreases like $`1/\mathrm{\Lambda }^2`$ if the 3-body coupling constant has the form $`g_3(\mathrm{\Lambda })=4mg_2(\mathrm{\Lambda })^2H(\mathrm{\Lambda })/\mathrm{\Lambda }^2`$, where the dimensionless function $`H(\mathrm{\Lambda })`$ has the form

$`H(\mathrm{\Lambda })={\displaystyle \frac{\mathrm{sin}[s_0\mathrm{ln}(\mathrm{\Lambda }/\mathrm{\Lambda }_{})\mathrm{arctan}(1/s_0)]}{\mathrm{sin}[s_0\mathrm{ln}(\mathrm{\Lambda }/\mathrm{\Lambda }_{})+\mathrm{arctan}(1/s_0)]}}`$ (3)
for some value of $`\mathrm{\Lambda }_{}`$. With this renormalization, 3-body observables have well-defined limits as $`\mathrm{\Lambda }\mathrm{}`$, but they depend on the parameter $`\mathrm{\Lambda }_{}`$ introduced by dimensional transmutation. Note that $`H(\mathrm{\Lambda })`$ is a periodic function of $`\mathrm{ln}(\mathrm{\Lambda })`$. Thus the renormalization of the field theory involves an ultraviolet limit cycle.

An effective field theory (EFT) of nucleons with contact interactions only has also been applied to the 3-nucleon problem Bedaque:1999ve ; Hammer:2001gh . This contact EFT involves an isospin doublet $`N`$ of Pauli fields with two independent 2-body contact interactions: $`N^{}\sigma _iN^cN^c\sigma _iN`$ and $`N^{}\tau _kN^cN^c\tau _kN`$, where $`N^c=\sigma _2\tau _2N^{}`$. Renormalization in the 2-body sector requires that the two coupling constants be adjusted as a function of the ultraviolet cutoff $`\mathrm{\Lambda }`$ to obtain the correct spin-singlet and spin-triplet scattering lengths $`a_s`$ and $`a_t`$. Renormalization in the 3-body sector requires the 3-body contact interaction $`N^{}\sigma _iN^cN^c\sigma _jNN^{}\sigma _i\sigma _jN`$ with a coefficient proportional to (3). Thus the renormalization involves an ultraviolet limit cycle.

The low-energy few-nucleon problem can also be described by an EFT that includes explicit pion fields as well as contact interactions between the nucleons. The renormalization of such an EFT does not involve any RG limit cycle. This suggests that the RG limit cycle in the EFT of Ref. Bedaque:1999ve  is just an artifact of the model. We will argue to the contrary that it is actually a hint of an infrared limit cycle in QCD.

Our argument is based on recent work in which an EFT with explicit pions was used to extrapolate nuclear forces to the chiral limit of QCD Beane:2002xf ; Epelbaum:2002gb . In this limit, the masses $`m_u`$ and $`m_d`$ of the up and down quarks are zero and the pion is an exactly massless Goldstone boson associated with spontaneous breaking of the chiral symmetry of QCD. According to the chiral extrapolation of Epelbaum, Meißner, and Glöckle Epelbaum:2002gb , the small binding energy 2.2 MeV of the deuteron is a fortuitous consequence of the physical values of $`m_u`$ and $`m_d`$. When extrapolated to the chiral limit, the deuteron has a much larger binding energy comparable to the scale $`1/(m_Nr_0^2)10`$ MeV set by the $`NN`$ effective range. Conversely, if extrapolated farther from the chiral limit, the deuteron’s binding energy decreases to 0 and then it becomes unbound.

This effect is illustrated in Fig. 1, which shows the chiral extrapolation of the inverse scattering lengths $`1/a_t`$ and $`1/a_s`$ as functions of $`m_\pi `$ from Ref. Epelbaum:2002gb . In the EFT of Ref. Epelbaum:2002gb , the coefficients of some of the 2-nucleon contact interactions are not well-constrained by the low-energy data. The bands in Fig. 1 are obtained by varying those coefficients over natural ranges. The width of the error band of course shrinks to zero at the physical value of $`m_\pi `$. The prediction of Ref. Epelbaum:2002gb  is that the critical value $`m_{\pi ,t}`$ at which $`1/a_t=0`$ is in the range 170 MeV $`<m_{\pi ,t}<`$ 210 MeV, which is not much larger than the physical value of $`m_\pi `$. The extrapolation of $`1/a_s`$ has larger uncertainties. It may increase to zero at some critical value $`m_{\pi ,s}`$ greater than 150 MeV, in which case the spin-singlet deuteron is bound for $`m_\pi >m_{\pi ,s}`$, or $`1/a_s`$ may remain negative. Beane and Savage Beane:2002xf  have argued that the extrapolation errors in the chiral limit are larger than estimated in Ref. Epelbaum:2002gb . For the small extrapolation to the region of larger $`m_\pi `$ where the deuteron becomes unbound, however, there is no controversy and the extrapolation errors are under control.

We now consider the chiral extrapolation of the 3-body spectrum. This could be calculated using an EFT with explicit pions. Alternatively, the chiral extrapolation can be calculated using the contact EFT of Ref. Bedaque:1999ve . The inputs required are $`a_s`$, $`a_t`$, and $`\mathrm{\Lambda }_{}`$ as functions of $`m_\pi `$, which can be calculated using an EFT with pions. For the inverse scattering lengths $`1/a_s`$ and $`1/a_t`$, we take the central values of the error bands obtained from the chiral extrapolation in Ref. Epelbaum:2002gb . The dependence of $`\mathrm{\Lambda }_{}`$ on $`m_\pi `$ could be determined from the chiral extrapolation of the triton binding energy using an EFT with pions, but this has not yet been calculated. Like the inverse scattering lengths, $`\mathrm{\Lambda }_{}`$ should vary smoothly with $`m_\pi `$. For small extrapolations of $`m_\pi `$ from its physical value, we can approximate $`\mathrm{\Lambda }_{}`$ by a constant. We use the value $`\mathrm{\Lambda }_{}=189`$ MeV for $`m_\pi =138`$ MeV obtained by taking the triton binding energy as the 3-body input.

In Fig. 2, we show the 3-body spectrum in the triton channel as a function of $`m_\pi `$. The crosses give the binding momenta $`\kappa =(mB_3)^{1/2}`$ of the physical deuteron and triton, while the dashed lines give the thresholds for decay into a nucleon plus a deuteron (left curve) or a spin-singlet deuteron (right curve) in the large-$`a`$ approximation. The circles indicate the triton ground state and excited state. In the region near $`m_\pi 175`$ MeV where the decay threshold comes closest to the 3-nucleon threshold $`\kappa =0`$, the excited state of the triton appears in the spectrum. The existence of this very shallow excited state is a hint that the system is close to an infrared limit cycle.

If, as in the case illustrated by Fig. 2, the critical values $`m_{\pi ,t}`$ and $`m_{\pi ,s}`$ at which $`1/a_t`$ and $`1/a_s`$ go to zero satisfy $`m_{\pi ,s}<m_{\pi ,t}`$, then the deuteron and spin-singlet deuteron are both bound in the region $`m_{\pi ,s}<m_\pi <m_{\pi ,t}`$. Since the decay threshold never extends all the way down to the 3-nucleon threshold $`\kappa =0`$, there cannot be an infrared limit cycle. We argue that it is possible, however, to tune the system to the critical RG trajectory for an infrared limit cycle by adjusting the up and down quark masses $`m_u`$ and $`m_d`$. In the next-to-leading order chiral extrapolation of Ref. Epelbaum:2002gb , only quark mass dependent operators proportional to $`m_u+m_d`$ enter. The extrapolation in $`m_\pi `$ can be interpreted as an extrapolation in the sum $`m_u+m_d`$, with $`m_um_d`$ held fixed. Changing $`m_um_d`$ changes the degree of isospin-symmetry breaking. Since the spin-singlet and spin-triplet channels correspond to isospin-triplet and isospin-singlet, respectively, the scattering lengths $`a_s`$ and $`a_t`$ will respond differently to changes in $`m_um_d`$. By tuning both $`m_u`$ and $`m_d`$, it should be possible to make $`1/a_t`$ and $`1/a_s`$ vanish simultaneously: $`m_{\pi ,t}=m_{\pi ,s}`$. This point corresponds to a critical RG trajectory for an infrared limit cycle. At this critical point, the triton has infinitely many increasingly-shallow excitations with an accumulation point at the 3-nucleon threshold. As the excitations get more and more shallow, the ratio of the binding energies of successive states rapidly approaches $`\lambda _0^2515`$. Only for the deepest states like the triton will there be significant deviations. Now consider a renormalization group transformation that integrates out energies above some scale $`\mathrm{\Lambda }^2/m_N`$. As $`\mathrm{\Lambda }`$ is decreased, the deepest 3-nucleon bound states are removed from the spectrum, leaving only those for which the deviations from the asymptotic ratio $`\lambda _0^2`$ are negligible. Thus a limit cycle with a discrete-scaling-symmetry factor $`\lambda _0`$ is approached in the infrared limit $`\mathrm{\Lambda }0`$.

The error bands in Ref. Epelbaum:2002gb  do not exclude the critical values $`m_{\pi ,t}`$ and $`m_{\pi ,s}`$ from satisfying $`m_{\pi ,t}<m_{\pi ,s}`$, in which case the deuteron and spin-singlet deuteron are both unbound in the region $`m_{\pi ,t}<m_\pi <m_{\pi ,s}`$. In this region, the decay threshold extends all the way down to the 3-nucleon threshold $`\kappa =0`$. At either of the two critical values $`m_{\pi ,t}`$ and $`m_{\pi ,s}`$, there is a 2-nucleon bound state at threshold. However there is no Efimov effect, because the 3-body problem is not sufficiently singular at short distances. The exact infrared limit cycle can only be obtained by using an additional tuning parameter such as $`m_um_d`$ to set $`m_{\pi ,t}=m_{\pi ,s}`$.

The infrared RG limit cycle of QCD has important implications for attempts to derive nuclear physics from lattice gauge theory. The computational effort for lattice simulations increases dramatically as the pion mass decreases and is prohibitive at the physical value. Lattice simulations are typically carried out at a value of $`m_\pi `$ that is 2-3 times larger than the physical value, and then a chiral extrapolation is made to to $`m_\pi =138`$ MeV. Unfortunately, this requires extrapolating past the region of $`m_\pi `$ where there is an RG limit cycle, which may introduce large extrapolation errors. The proximity of physical QCD to the critical trajectory for the infrared limit cycle explains the success of Efimov’s program Efi81  for describing the 3-nucleon problem using zero-range forces or a contact EFT Bedaque:1999ve . The apparent convergence of the effective-range corrections for momenta of order $`m_\pi `$ Hammer:2001gh  can be explained if the momentum expansion in the contact EFT is in powers of $`p/m_\pi ^{}`$, where $`m_\pi ^{}175`$ MeV is the critical value of the pion mass, which is significantly larger than the physical value. Our final conclusion is that the ultraviolet limit cycle in the contact EFT Bedaque:1999ve  is not just an artifact of the model, but a hint of the existence of the infrared limit cycle of QCD.

This research was supported in part by DOE grant DE-FG02-91-ER4069. We thank U.-G. Meißner for useful discussions, E. Epelbaum for providing the scattering lengths $`a_s`$ and $`a_t`$ as functions of $`m_\pi `$, and R.J. Furnstahl for comments on the manuscript.