text:
# Water Vapor in Carbon-rich AGB Stars from the Vaporization of Icy Orbiting Bodies

## 1 Introduction

One of the most exciting developments in astronomy during the past decade has been the unequivocal detection of extra-solar system planets of size as small as a few Jupiter masses (e.g. Mayor & Queloz 1995, Marcy & Butler 1996). This development suggests the possibility that planetary systems around other stars harbor yet smaller bodies, including Earth-sized planets, asteroids and comets. The presence of comets, in particular, has been suggested by the $`\beta `$ Pictoris phenomenon, in which time-variable line absorption is observed in a stellar spectrum at substantial Doppler shifts relative to the central star (Ferlet, Vidal-Madjar, & Hobbs 1987), suggesting comets evaporating as they speed near the central star. The discovery of extra-solar planetary systems also raises the question of how such systems change as the stars that they surround evolve.

Stern, Shull & Brandt (1990; hereafter SSB90) have pointed out when a star evolves off the main sequence, the resultant increase in its luminosity must lead to the vaporization of any icy bodies orbiting within several hundred AU, an effect that could result in the release of significant amounts of water vapor if the star were surrounded by a Kuiper Belt or Oort cloud. SSB90 considered the vaporization of icy bodies of various radii at various distances from a post-main-sequence star of constant assumed luminosity $`6000\mathrm{L}_{}`$, and computed the resultant rate at which water vapor is deposited in the circumstellar environment. They suggested that this process could be responsible for the large water vapor and OH abundances typically observed in the circumstellar envelopes of oxygen-rich stars.

Although large water vapor abundances are in any case expected around oxygen-rich stars, and although the water outflow rates around many such sources greatly exceed what can plausibly be explained by the SSB90 model, recent measurements of modest water vapor abundances in two carbon-rich circumstellar outflows (Herpin & Cernicharo 2000, Melnick et al. 2001) motivate renewed attention to the mechanism proposed by SSB90. In particular, Melnick et al. have used the Submillimeter Wave Astronomy Satellite (SWAS) to detect water vapor emission from the classic carbon star IRC+10216. Carbon stars, which are late-type AGB stars where the carbon-to-oxygen ratio is greater than $`1`$, are expected to harbor almost no water vapor in their circumstellar environments; the equilibrium chemistry of oxygen is dominated by CO and there is ordinarily very little oxygen left to form any other molecules. The average water abundance predicted by standard chemical models for IRC+10216 (Millar et al. 2000) is less than $`10^{12}`$ (Millar & Herbst 2001) and thus the presence of detectable water vapor in IRC+10216 may be a distinctive signature of the vaporization of orbiting icy bodies.

The existence of this signature rests critically upon the survival of such bodies. In particular, it is not possible to determine from the study of SSB90 whether any icy bodies will survive the $`1`$ Gyr of post-main-sequence evolution that precedes the carbon-rich phase (at least at radii where such bodies are plausibly present in the first place and where they would be subject to vaporization by a carbon-rich AGB star). To address this question we have extended the work of SSB90 by considering (1) the exact luminosity variations expected in a post-main-sequence star; (2) the size distribution and total mass of icy bodies expected if extra-solar cometary systems are similar to our own Kuiper Belt (the properties of which have been greatly elucidated \[e.g. Jewitt & Luu 2000\] since the work of SSB90); and (3) the specific application to carbon-rich stars, where the vaporization of icy bodies has the most distinctive observational signature. Our paper is organized as follows: in §2, we detail our theoretical model of the vaporization of a Kuiper Belt analog around a post-main-sequence star; in §3 we present the results of our calculations; the results are discussed in §4.

## 2 Calculations

For our standard model, we considered a star of mass $`M_{}=1.5M_{}`$ and followed its luminosity evolution from the main sequence to the late AGB stage using evolutionary tracks kindly provided by Allen Sweigart. The luminosity, $`L`$, during the last $`2`$Myr of the AGB track (the thermal pulsation or “TP-AGB” phase) is shown in Figure 1; the zero of time corresponds to the first thermal pulse on the AGB. The stellar evolution model assumes a mass-loss rate, $`\dot{M}_{}`$, given by the expression of Reimers (1975), $`(\dot{M}_{}/M_{}\mathrm{yr}^1)=4\times 10^{13}\eta _R(L/L_{})(R_{}/R_{})(M_{}/M_{})^\mathit{1}`$, where $`R_{}`$ is the stellar radius, and with the “Reimers parameter”$`\eta _R`$ equal to 0.4. With this prescription for the mass-loss rate, the stellar mass has decreased to $`1.4M_{}`$ by the start of the AGB phase and to $`0.92M_{}`$ by the end of the AGB phase.

Using this model for the luminosity evolution of a $`1.5M_{}`$ star, we have calculated how a population of orbiting icy bodies would evolve. The mass loss rate per unit surface area, $`\dot{m}(T)`$, is described by the equation $`\dot{m}(T)=P_s(T)(\mu /2\pi kT)^{1/2}`$, where $`\mu `$ is the molecular mass of water, $`k`$ is the Boltzmann constant and $`P_s(T)`$ is the vapor pressure of water ice at temperature $`T`$. This equation, which neglects gravitational effects, applies to bodies of radius $`\stackrel{<}{}1000`$ km for which the escape velocity is smaller than the typical outflow velocity of the vaporizing molecules ($`1\mathrm{km}\mathrm{s}^1`$). Furthermore, it applies strictly to bodies composed of pure water ice rather than to bodies composed of a mixture of ice and refractory material. In large icy bodies (the size of Pluto and Charon), it is very plausible that differentiation will lead to a structure in which relatively pure ice surrounds a rocky core (McKinnon, Simonelli & Schubert 1997). In smaller bodies, however, the dust and ice will be well mixed and the effects of the refractory material on the vaporization rate is less clear – and may depend critically upon the size distribution of the dust (Prialnik 1992; Orosei et al. 1995). The average surface temperature in equilibrium, $`T`$, reached by an icy body at distance $`R`$ from a star of luminosity $`L`$ is given by $`ϵ\sigma T^4=[(1A)L/(16\pi R^2)H\dot{m}(T)]`$, where $`A`$ is the albedo, $`H=2.45\times 10^{10}\mathrm{erg}\mathrm{g}^1`$ is the specific heat of vaporization (Prialnik 1992), $`ϵ`$ is the emissivity and $`\sigma `$ is the Stefan-Boltzmann constant. Following SSB90, we adopt an albedo $`A=0.04`$, an emissivity $`ϵ=1`$, and a vapor pressure $`P_s(T)=10^{9.1838372403.4/T}\mathrm{torr}`$ (Lebofsky 1975).

The mass-loss caused by vaporization reduces the radius of each orbiting icy body. For a spherical body of uniform density $`\rho `$, the change in radius after time $`t`$ is given by $`\mathrm{\Delta }r(t)=_0^t(\dot{m}(t^{})/\rho )𝑑t^{}`$. After time $`t`$, all objects of initial radius smaller than $`\mathrm{\Delta }r`$ have been completely vaporized. In computing the total mass loss rate from a population of orbiting icy bodies, we assume an initial size distribution of the form $`dn/drr^q`$, where $`dn`$ is the number of objects of radius $`r`$ to $`r+dr`$. Recent observations (Jewitt & Luu 2000) of our own Kuiper Belt suggest $`q4`$ for radii $`r`$ between $`r_{\mathrm{min}}1`$km and $`r_{\mathrm{max}}1200`$km (the size of Pluto); in other words, the size distribution has equal amounts of ice mass per logarithmic radius interval and the available reservoir of ice is dominated neither by small nor by large bodies. Although well-motivated by observations of the past decade, our assumption that the Kuiper belt contains objects with a range of sizes is not a critical one, and calculations with a single assumed object size yield qualitatively similar results.

## 3 Results

Our calculation shows that at the astrocentric radii of the classical Kuiper Belt, $`R=3050\mathrm{AU}`$ (Jewitt & Luu 2000), even large icy bodies are entirely vaporized prior to the TP-AGB phase111Note, however, that any icy bodies originally at the radius of the classical Kuiper Belt may have experienced significant outwards migration during the post-main-sequence evolution as a result of stellar mass-loss.. Thus by the time the dredge-up of carbon has made the photospheric C/O ratio greater than 1, the release of water vapor can only take place if icy bodies are present at larger radii. Figure 2 shows the “complete vaporization radius”, $`R_{\mathrm{vap}}`$, within which every icy body has been destroyed, as a function of time after the start of the TP-AGB phase. The solid curve shows the results for a size distribution extending to $`r_{\mathrm{max}}=1200`$km, the dotted curve for the case $`r_{\mathrm{max}}=120`$km, and the dashed curve for the case $`r_{\mathrm{max}}=12`$km. Clearly, the complete vaporization radius is a slowly increasing function of $`r_{\mathrm{max}}`$, because larger bodies have a smaller surface-area-to-mass ratio and thus a longer lifetime.

Figure 3 shows the evolution of the water outflow rate and abundance for a collection of bodies located in circular orbits at fixed distance $`200`$AU from the central star. The mass-loss rate per unit initial ice mass, $`\dot{M}(\mathrm{H}_2\mathrm{O})/M_0(\mathrm{ice})`$, is shown in Figure 3a for the case $`r_{\mathrm{max}}=1200`$km. The quantity $`M_0(\mathrm{ice})`$ refers to the total mass of water ice present in the collection of orbiting bodies prior to the onset of evaporation. The same results appear in Figure 3b, but are now expressed as a water abundance, $`x(\mathrm{H}_2\mathrm{O})n(\mathrm{H}_2\mathrm{O})/n(\mathrm{H}_2)`$. Given a water mass outflow rate of $`\dot{M}(\mathrm{H}_2\mathrm{O})`$, the water abundance in the outflowing gas is $`\frac{1}{9}\times \dot{M}(\mathrm{H}_2\mathrm{O})/\frac{3}{4}\times \dot{M}_{}`$ relative to H2, the factor of $`\frac{3}{4}`$ being the mass fraction of hydrogen and the factor of $`\frac{1}{9}`$ being the ratio of the H2 molecular mass to that of H2O. Because the H2 mass-loss rate is proportional to the Reimers parameter, $`\eta _R`$, the quantity plotted in Figure 3b is $`x(\mathrm{H}_2\mathrm{O})\eta _R/M_0(\mathrm{ice})`$. We note that our use of the Reimers parameter $`\eta _R`$ in Figure 3b is merely a useful way of parameterizing the mass-loss rate and that the results presented here apply regardless of whether “Reimers Law” is correct – i.e. regardless of whether $`\eta _R`$ is a constant. Indeed, there is a developing consensus (e.g. Blöcker 1995, Willson 2000) that $`\eta _R`$ is an increasing function of luminosity and that mass-loss occurs preferentially near the tip of the AGB. The current value of $`\eta _R`$ for the high mass-loss star IRC+10216 is $`10`$.

In Figures 4a and 4b, the water outflow rate and abundance are now shown as a function of astrocentric distance, $`R`$, at the three example evolutionary stages marked by open symbols on Figures 1 and 3. As in Figure 3, these results apply to a collection of icy bodies all located at a single distance, $`R`$, from the star. At any given time, the water outflow rate per unit initial ice mass shows a peak at some particular distance. At small distances from the star, the outflow rate and abundance increase as the assumed distance gets larger because the amount of icy material remaining is an increasing function of $`R`$. (This increase is a gradual one because we assume the presence of a range of sizes for the orbiting icy bodies.) At large astrocentric distances, the outflow rate and abundances drop with increasing $`R`$ because of the diminishing stellar flux. The radius at which $`\dot{M}(\mathrm{H}_2\mathrm{O})/M_0(\mathrm{ice})`$ peaks moves outwards as the star evolves because of the steady increase in stellar luminosity on the AGB. Figure 4a implies that for icy bodies at the right distance from the star, water outflow rates as large as $`10^5M_0(\mathrm{ice})`$ yr-1 can plausibly be attained during the TP-AGB phase.

## 4 Discussion

The results presented above can now be discussed in the context of recent SWAS observations of the carbon-rich AGB star IRC+10216. The extraordinary stability of the SWAS receivers has allowed radiometric performance to be achieved in observations of duration as long as 200 hours (Melnick et al. 2000). Long duration observations have thereby achieved sensitivities considerably better than had been envisaged prior to the launch of SWAS. Such observations have led to the unexpected detection of the $`1_{10}1_{01}`$ 557 GHz transition of water vapor in emission toward IRC+10216 (Melnick et al. 2001; hereafter M01). The measured line strength implies a water abundance in the outflowing circumstellar envelope in the range 4 – 24 $`\times 10^7`$ (M01) corresponding to a water mass-loss rate of $`24\times 10^{10}M_{}`$ yr-1 $`=0.61.4\times 10^4M_{}`$ yr-1. These exact values depend upon the H2 mass-loss rate and the distance to IRC+10216, both of which are somewhat uncertain.

The exact initial mass of IRC+10216 is poorly-known but is believed to be between $`1.5`$ and $`4M_{}`$ (Forestini, Guelin & Cernicharo 1997, Kahane et al. 2000). We have carried out calculations analogous to those presented in §3 with a stellar mass of $`4M_{}`$ in place of the standard case $`M_{}=1.5M_{}`$; the results are qualitatively similar to those for the standard case, implying that our results are neither critically dependent upon the mass of the star nor upon the exact (and somewhat uncertain) details of the luminosity evolution model.

The results shown in Figure 4a and discussed in §3 above imply that the vaporization of icy bodies can explain the observed abundances of water vapor in the IRC+10216 outflow, given a total initial ice mass, $`M_0(\mathrm{ice})`$, that can be as little as $`10M_{}`$ if the icy bodies are located at a suitable radius. This value increases to $`50M_{}`$ in a distributed model in which the icy bodies are assumed to be spread over a broad range of radii (40 – 400 AU) with the surface density of icy bodies proportional to $`R^2`$. The exact value of $`M_0(\mathrm{ice})`$ required depends upon several uncertain parameters, including the actual water abundance, the range of astrocentric radii at which icy bodies are present, the maximum object radius, and the evolutionary stage222From the fact that IRC+10216 possesses one of the highest mass-loss rates of any known carbon star – a mass-loss rate some two orders of magnitude above the median value for carbon AGB stars – it seems likely that this source is in a stage of enhanced luminosity close to the tip of the AGB. of the star.

Although the properties of any Kuiper Belt analog around IRC+10216 might be quite different from those of our own Kuiper Belt, it is of interest to compare the required value of $`M_0(\mathrm{ice})`$ with the mass inferred for the one system of which we have direct knowledge. Dynamical considerations place a limit of $`1M_{}`$ upon the mass located at heliocentric radii in the range 30 – 50 AU (i.e. in the ”classical Kuiper belt”), while recent detections of large Kuiper Belt Objects (KBOs) suggest a total mass $`0.1M_{}`$ (Jewitt 1999), of which $`50\%`$ is likely water ice. Thus the observed amount of water ice currently present in the nearby classical Kuiper belt is far below what would be needed to account for the water abundance observed in IRC+10216. However, models for the formation of large KBOs such as Pluto (Kenyon & Luu 1998), suggest that the Kuiper Belt mass must $`\mathrm{𝑜𝑟𝑖𝑔𝑖𝑛𝑎𝑙𝑙𝑦}`$ have been at least $`10M_{}`$ during the early history of the solar system; thus the vast majority of the material has likely been ejected to larger heliocentric distances in the Kuiper belt or Oort cloud where it cannot be detected in current surveys.

In addition to water vapor, the vaporization of comet-like icy bodies could release a distinctive mixture of other oxygen-bearing molecules that are not expected in carbon-rich environments but which are typically observed in comets. These include carbon dioxide and methanol. We are not aware of any detections of these molecules333Latter & Charnley (1996) reported the detection of millimeter line emission at the frequencies of several methanol transitions, but subsequently favored C4H and C4H2 as the species responsible for the emission originally attributed to methanol. toward IRC+10216 nor of any upper limits that provide useful constraints upon their abundances.

The discussion presented above suggests that the vaporization of orbiting icy bodies is a plausible explanation for the water vapor detected in the circumstellar environment of IRC+10216. In addition, the vaporization scenario makes a specific prediction for the spatial distribution of the circumstellar water vapor. As the results in Figure 2 show, the vaporization of icy bodies injects water vapor into the outflow only at radii greater than $`75`$ AU. As noted by M01, the strength of the low-lying submillimeter transition observed by SWAS relative to higher-excitation far-infrared transitions will be higher in the vaporization scenario than would be the case were the water vapor abundance uniform throughout the envelope. While upper limits on far-infrared water line luminosities obtained by the Infrared Space Observatory (Cernicharo et al. 1996) provide marginal evidence against a uniform water abundance (M01), definitive observations of the water line ratios will have to await the greater sensitivity of the Herschel Space Observatory (HSO, formerly known as FIRST). The expected sensitivity of HSO – probably three orders of magnitude better than SWAS for observations at 557 GHz – will also allow searches for water vapor to be carried out in a much larger sample of carbon-rich AGB stars.

We thank Allen Sweigart for providing us with evolutionary tracks of 1.5 and 4 $`M_{}`$ stars. We gratefully acknowledge helpful discussions with Gary Melnick, Eric Herbst, Al Glassgold and David Hollenbach. This work was supported by subcontract SV252005 from the Smithsonian Institution.

Figure Captions

Fig. 1 – Evolutionary tracks for a 1.5 $`M_{}`$ star (kindly provided by Allen Sweigart). The luminosity, $`L`$, during the last $`2`$Myr of the AGB track (the thermal pulsation or “TP-AGB” phase) is shown in Figure 1; the zero of time corresponds to the first thermal pulse on the AGB. Open symbols mark the three example evolutionary stages plotted in Figure 4.

Fig. 2 – Vaporization radius, $`R_{\mathrm{vap}}`$, within which every icy body has been destroyed, as a function of time after the start of the TP-AGB phase. The solid curve shows the results for a size distribution extending to $`r_{\mathrm{max}}=1200`$km, the dotted curve for the case $`r_{\mathrm{max}}=120`$km, and the dashed curve for the case $`r_{\mathrm{max}}=12`$km

Fig 3. – Evolution of the water outflow rate and abundance for a collection of bodies located in circular orbits at fixed distance $`200`$AU from the central star. The mass-loss rate per unit initial ice mass, $`\dot{M}(\mathrm{H}_2\mathrm{O})/M_0(\mathrm{ice})`$, is shown in Figure 3a for the case $`r_{\mathrm{max}}=1200`$km. The quantity $`M_0(\mathrm{ice})`$ refers to the total mass of water ice initially present in the collection of orbiting bodies. The same results appear in Figure 3b, but are now expressed as a water abundance, $`x(\mathrm{H}_2\mathrm{O})n(\mathrm{H}_2\mathrm{O})/n(\mathrm{H}_2)`$. Open symbols mark the three example evolutionary stages plotted in Figure 4.

Fig 4. – Water outflow rate and abundance as a function of astrocentric distance, $`R`$, at the three example evolutionary stages marked by open symbols on Figures 1 and 3. The short-dashed line refers to the evolutionary stage marked by a triangle in Figures 1 and 3, the long-dashed line to that marked by a diamond, and the solid line to that marked by a star. As in Figure 3, the results in Figures 4a and 4b apply to a collection of icy bodies all located at a particular distance, $`R`$, from the star.


format: text

text:
# Fast Clustering Analysis of Inhomogenous Megapixel CMB Maps

## 1 Introduction

With the successful launch of MAP and the advancing Planck schedule, megapixel Cosmic Microwave Background (CMB) maps are just around the corner. To fully realize their potential of constraining cosmological parameters within a few percent (e.g., Spergel 1994; Knox 1995; Hinshaw, Bennett, & Kogut 1995; Jungman et al. 1996; Zaldarriaga, Spergel & Seljak 1997; Bond, Efstathiou, & Tegmark 1997), data analysis methods have to face hitherto unprecedented challenges. The standard maximum likelihood developed for COBE (e.g., Górski 1994; Górski et al. 1994, 1996; Bond 1995; Tegmark & Bunn 1995; Tegmark 1996; Bunn & White 1996; Bond, Jaffe, & Knox 1998, 2000) is already pushing supercomputers for balloon-borne and ground based experiments with $`N10^5`$ pixels (de Bernardis et al. 2000; Netterfield et al. 2001; Hanany et al. 2000; Jaffe et al. 2000; Martin et al. 1996; Miller et al. 1999; Peterson et al. 2000), and it is clearly inadequate for analysing megapixel maps on any future supercomputer. This letter addresses an important step in the full data processing pipeline from component separation/mapmaking to cosmological parameter estimation: the fast estimation of the angular power spectrum, $`C_{\mathrm{}}`$ from megapixel map with realistic geometry (sky cut, cut out holes), and inhomogeneous coverage and noise.

A recent surge of activity motivated by the obvious need for fast alternatives to the standard maximum likelihood estimation produced an array of promising techniques. Several of them use symmetries of a particular experiment to gain speed. An experiment specific technique for MAP was developed by Oh, Spergel and Hinshaw (1998) making use of the planned approximate azimuthal symmetry of the coverage/noise. Similar developments are under way for Planck (e.g., Wandelt 2000) possibly exploiting symmetries in its scanning strategy. Other techniques give up optimality in favor of speed. Correlation functions with simple weights were advocated by Szapudi et al. (2001, hereafter SPPSB), while an analogous method based on direct spherical Fourier transform was developed by Hivon et al. (2001). In the standard CMB formalism (e.g. Bond et al. 1998) both of these can be thought of as quadratic estimators with suboptimal weights, yet, they tend to produce results very close to optimal.

This letter presents significant improvements of the correlation function method. We present a fast algorithm for extracting the pixel space estimator using fast spherical harmonics transform, introduce noise weighting, test the use of Monte Carlo realizations to calculate noise correlation functions in the estimator, and finally illustrate the possibility of additional signal weighting.

## 2 Summary of the Method

Our latest recipe to extract $`C_{\mathrm{}}`$’s from a large CMB map is an improved version of SPPSB: extract the two-point correlation function with an unbiased weighted estimator sampled at the roots of Legendre polynomials, then integrate with a Gauss-Legendre quadrature to obtain the $`C_{\mathrm{}}`$’s. Next we review the correlation function estimator.

Let us denote the temperature fluctuations at a sky vector $`n`$, a unit vector pointing to a pixel on the sky, with $`T(n)`$. In isotropic universes the two-point correlation function is a function only of the angle between the two vectors and can expanded into a Legendre series,

$$\xi _{12}=T(n_1)T(n_2)=\underset{\mathrm{}}{}\frac{\mathrm{}+1/2}{\mathrm{}(\mathrm{}+1)}𝒞_{\mathrm{}}P_{\mathrm{}}(\mathrm{cos}\theta ),$$
(1)
where $`\mathrm{cos}\theta =n_1.n_2`$ is the dot product of the two unit vectors, $`P_{\mathrm{}}(x)`$ is the $`\mathrm{}`$-th Legendre polynomial, and the $`𝒞_{\mathrm{}}`$ coefficients realize the angular power spectrum of fluctuations. If the CMB anisotropy is Gaussian, which is expected to be an excellent approximation, the correlation function, or the $`𝒞_{\mathrm{}}`$’s, yields full statistical description.

In reality each pixel value, $`\mathrm{\Delta }_i`$, contains contributions from the CMB and noise; the latter is also assumed to be Gaussian with a correlation matrix (the noise matrix) $`N_{ij}`$, determined during map-making. The full pixel-pixel correlation matrix is $`C_{ij}=\xi _{ij}+N_{ij}.`$ matrix, as a sum of noise correlation functions measured in a set of noise realizations. Therefore we use the modified version of the estimator by SPPSB

$$\stackrel{~}{\xi }(\mathrm{cos}\theta )=\underset{ij}{}f_{ij}(\mathrm{\Delta }_i\mathrm{\Delta }_j\frac{1}{M}\underset{k=1}{\overset{M}{}}n_i^kn_j^k),$$
(2)
where $`n_i^k`$ is one of $`M`$ realizations of the noise for pixel $`i`$. Our estimator can be extracted efficiently and accurately for a general noise matrix, as long as a fast noise generator exists. The weights $`f_{ij}`$ merit further discussion. For $`f_{ij}(S+N)^2`$ the above corresponds to an optimal quadratic estimator (where $`S`$ denotes the matrix corresponding to $`\xi _{ij}`$). Unfortunately, the cost to calculate the optimal weights would be prohibitive as it is for any other version of the optimal quadratic estimator; speed is gained by using heuristic weights instead of the optimal ones. SPPSB used window weighing, $`f_{ij}=0`$ unless the pair of pixels belong to a particular bin in $`\mathrm{cos}\theta `$, and $`_{ij}f_{ij}=1`$. The analysis of the next section uses more general noise weighting. Other possibilities exist and are worth to explore, however, they can only affect slightly the error bars of our unbiased estimator. The next section demonstrates that heuristic noise weighting produces results close to the theoretical expectations. Heuristic window weighting is natural in pixel space where the window is diagonal, and it amounts to edge effect correction (Szapudi & Szalay 1998). Heuristic noise weighting takes this idea one step further, and it is again natural in pixel space as long as the noise is localized as it is for MAP. It is worth reemphasizing that azimuthal symmetry is not required.

## 3 Fast Estimation of the Correlation Function at Legendre Roots

The novel correlation function estimator is based on fast spherical harmonics decomposition; our aim is to make use of their speed, based on the FFT’s of isolatitude pixel annuli. Practical implementation of our estimator is based on the fast transforms available in the HEALPix 555http://www.eso.org/kgorski/healpix package. Let us define two unit vectors with angle $`\theta `$ apart, e.g., $`n_1=(0,0)`$, and $`n_2=(\theta ,0)`$ in spherical coordinates. In the continuum limit, the raw weighted pairwise estimator is an average over the rotation group

$$D_{f\delta }(cos\theta )=\frac{1}{8\pi ^2}𝑑gf\delta (R(g)n_1)f\delta (R(g)n_2),$$
(3)
where the integration is over SO(3), $`8\pi ^2`$ is the full volume of the Lie group, and R(g) is the transformation corresponding to the Lie group element. The sum in equation (2) is replaced by an integral, multiplicative weighting is assumed: $`f_{ij}=f_if_j`$, and $`f\delta f_i\mathrm{\Delta }_i`$ (both taken at the same vector), corresponds to the weighted temperature fluctuations in the continuum limit. After a simple calculation based on the general orthogonality of group representations

$$D_{f\delta }(cos\theta )=\underset{lm}{}a_{lm}^2\frac{1}{4\pi }P_{\mathrm{}}(cos\theta ),$$
(4)
where the $`a_{lm}`$’s are the coefficients of the spherical harmonics decomposition of $`f\delta `$. The ensemble average of the equation yields $`D_{f\delta }=(\xi +N)D_f`$, where we introduced the weight correlation function $`D_f`$, calculated similarly to $`D_{f\delta }`$. The unbiased correlation function estimator can be obtained as

$$\stackrel{~}{\xi }(cos\theta )=\frac{D_{f\delta }(cos\theta )\overline{D}_{fn}(cos\theta )}{D_f(cos\theta )},$$
(5)
where $`\overline{D}_{fn}`$ is the average raw noise correlation function calculated from $`M`$ simulations; obviously $`\stackrel{~}{\xi }=\xi `$. According to equation (4) this estimator is calculated via two spherical harmonics transforms in discrete pixel space, with the summation performed up to $`\mathrm{}_{max}`$. This is a fast realization of the estimator in equation (2) with multiplicative weights.

We cross-compared the above algorithm with a naive $`N^2`$ correlation code, and the results were identical. To obtain the angular power spectrum, we then sample the above equation exactly at the roots of Legendre polynomials, and perform Gauss-Legendre integration, as done by SPPSB.

The expected scaling is $`N^{3/2}\mathrm{log}N`$, corresponding to an effective $`N^{1.61}`$ slope for the present implementation in the range of $`N`$ we tested. The speed between calculating the correlation estimator and the $`C_l`$ inversion is divided approximately evenly; on a 500Mhz Dec alpha our present HEALpix compatible implementation takes $``$260 seconds, Planck analysis is feasible in about 40 minutes.

The above procedure appears unnatural at first: if we are using FFT’s, and our final goal is to estimate the angular power spectrum, why do we leave Fourier space at all? The reason is simple: window and noise are localized in pixel space, thus we can construct an unbiased estimator with a simple normalization. Such a normalization in pixel space is equivalent to an edge effect correction, and heuristic weights can be constructed intuitively and naturally. In principle an analogous procedure can be designed in $`\mathrm{}`$ space, but the non-diagonality of the window function results in a complex coupling matrix infused with Wigner $`3j`$ symbols, as Hivon et al. (2001) have shown in a tour de force calculation. The computation and inversion of this matrix to unbias the estimator is highly non-trivial numerically, especially in the presence of noise and complex geometry; so far they have demonstrated numerical feasibility for a relatively simple ellipsoidal window without any noise weighting. It is not inconceivable that this procedure can be successfully extended to more general windows and noise weighting, but at a price which can be regarded as unnecessary complexity; therefore we recommend the simple technique of constructing an unbiased weighted pixel space estimator instead.

## 4 Application: MAP Simulations

We have generated 1200 MAP simulations using HEALPix with inhomogenous sinusoidal coverage assuming a detector sensitivity of $`20\mu K`$ for a $`0.3^{}\times 0.3^{}`$ pixel, taken from the MAP homepage, (http://map.gsfc.nasa.gov/). We assumed a Gaussian beam of $`18`$ arcminutes, we neglected any sidelobes or asymmetries. Our simulations are similar, although not perfectly equivalent to that of OSH. In addition to the noisy MAP simulations, we have generated 1200 pure noise simulations to be used in our estimator of Equation (2). Our noise weighting was inversely proportional to the expected noise in a pixel, $`f_{ij}1/(\sigma _{n_i}^p\sigma _{n_j}^p)`$, with $`p=1`$ motivated by prewhitening, and $`p=2`$ by approximate minimum variance estimator. Both performed quite similarly, we show the results from $`p=2`$. We used a galactic cut of $`\pm 20`$ degrees.

The noise correlations were calculated with the code and were subtracted from the correlation function according to equation (2). While the noise was assumed to be diagonal, our practical implementation does not make use of this property, except in the heuristic noise weighting we adopted. If this assumption would break, our method would not be affected, as long as the simulations can be generated quickly.

The $`C_{\mathrm{}}`$ integration was performed via a Gauss-Legendre quadrature. Beam and pixel window functions were accounted for as usual. The results are displayed in Figure 3. Along with the average $`C_{\mathrm{}}`$, the theoretical input, as well as theoretical and measured errors are displayed, the former from the following approximation (Hivon et al. 2001)

$$\mathrm{\Delta }𝒞_{\mathrm{}}=\sqrt{\frac{2}{\nu _l}}\left(𝒞_{\mathrm{}}+B(\mathrm{})^2W(\mathrm{})^2N(l)\right),$$
(6)
where $`\nu _l=(2\mathrm{}+1)f_{sky}w_2^2/w_4`$ with $`w_m=\frac{1}{4\pi }𝑑\mathrm{\Omega }f^m(\mathrm{\Omega })`$, the $`m`$-th moments of the weight function $`f_i`$, $`B(\mathrm{}),W(\mathrm{}),N(\mathrm{})`$ represent beam, pixel, and noise effects, respectively, and $`f_{sky}`$ is the sky fraction in the map. According to Figure 3 our method is unbiased, and accurate for MAP: the error bars track closely the theoretical expectation within the accuracy of the approximate theory. The suboptimality of our technique is expected to be small for MAP-like surveys, but this needs to be confirmed by detailed comparison to optimal techniques. Our figure is similar to Figure 4 of Oh et al. (1999), but there are small differences: it appears that we have assumed a slightly higher noise and/or wider beam. It would be worthwhile to perform direct comparison of the two methods to assess the degree of suboptimality of ours, and the sensitivity of the experiment specific technique to a potential breakdown of the underlying assumptions.

## 5 Summary and Discussions

We have demonstrated that i) correlation functions with generalized edge correction can be extracted with full accuracy for a megapixel map in minutes on a modest personal workstation ii) the resulting $`C_{\mathrm{}}`$ computed from Monte-Carlo realizations of MAP observations with inhomogeneous coverage agree quite well with the (approximate) theoretical error bars inferred from the moments of the window function applied to the data to construct our estimator, as well as the cosmic variance. The results should be very close to optimal at high $`\mathrm{}`$ where the noise variance become dominant provided that the noise matrix is diagonally dominant, but this needs to be confirmed by a detailed comparison with optimal methods. The theoretical scaling of our method is $`N^{3/2}\mathrm{log}N`$, with a measured scaling of $`N^{1.61}`$ in our present implementation, for the range of $`N`$ we tested.

In the previous tests we have used MAP-like diagonal noise/coverage, no other (e.g., azimuthal) symmetries were assumed about noise, or geometry of the map. Cut out holes around bright sources, or any irregularity in the sampling, represent only minor perturbations to our method, therefore they should have no discernible effect on its speed or performance; the investigation of this point in sufficient detail is left for future work. Our approach is straightforward to generalize for non-diagonal noise matrix, specially in the case where the off-diagonal elements of the noise matrix correspond primarily to pairs of pixels of fixed angular separation, as is expected for MAP due to the differential nature of the measurement. In the general case, our main assumption is the fast generation of noise realizations in pixel space, which will be supplied by fast mapmaking methods (e.g. Wright 1996; Prunet et al. 2001).

At the heart of our method is the use of heuristic weighting corresponding to general edge effect correction. (Indeed, noise can be considered as a fuzzy window, or window can be considered as infinite noise). The simplicity of our weighting amounts to a substantial gain in speed over the costly optimal weighting $`(S+N)^2`$. In addition, the use of spherical harmonics decomposition facilitates heuristic signal weighting, which is natural in $`\mathrm{}`$ domain. The sum of the $`a_{lm}`$’s in equation (4) can be recognized as a pseudo $`C_l`$ of the weighted fluctuations estimated with uniform weights in $`\mathrm{}`$ space. For low $`\mathrm{}`$’s, for complicated noise and window patterns, it might be desirable to sum the $`a_{lm}`$’s with weights obtained from the correlation matrix $`<a_{l_1m_1}a_{l_2m_2}>`$ evaluated from Monte Carlo simulations. This matrix can be diagonalized e.g. up to a certain $`\mathrm{}`$, or within a band of $`\mathrm{\Delta }\mathrm{}`$, to calculate approximate (semi-heuristic) weights for our estimator of the correlation function. While such improvements do not seem necessary for MAP, a clear upgrade path for our method exists for more complex experiments if needed. Although more natural in $`\mathrm{}`$ space, heuristic signal weighting is possible in pixel space as well (e.g., Colombi, Szapudi, & Szalay 1998).

Our technique has further potential besides speeding up $`C_l`$ estimation: it opens up a full array of possible applications and generalizations. These include cross-correlations (between channels for component separation, between LSS and CMB, B-type polarization and lensing, etc) non-Gaussianity (e.g. 3-point function/bispectrum, cumulant correlators for SZ, lensing, etc), vector and tensor correlations (for polarization). The idea of real space window and noise weighting is applicable to power spectrum estimators of galaxies, clusters, lensing etc. as well. A Euclidian version of our algorithm is entirely analogous to the spherical case. It will be useful for fast edge and noise corrected estimation of the power spectrum, bispectrum, $`N`$-point correlation function and cumulant correlators in galaxy catalogs. These generalizations presently under implementation will be discussed in subsequent papers and included in a later version of SpICE.

We would like to thank Dick Bond, Dmitry Pogosyan, Alex Szalay for their help, and Carlo Contaldi and Eric Hivon for useful discussions. This research was supported by NASA through the Applied Information Systems Research program (NAG5-10750).


text:
# Combinatorial identities for binary necklaces from exact ray-splitting trace formulae

## I Introduction

Wave propagation in systems with sharp interfaces is a fundamental problem in the natural sciences and engineering. Well-known examples include light waves impinging on a water-air interface or sound waves propagating in layered media . All these systems have one feature in common: the splitting of the incident wave into reflected and transmitted components. In the geometrical optics limit of small wave lengths, incident, reflected and transmitted waves are described by rays. The rays are split at the interface; hence the name “ray-splitting systems”  for the whole class of wave systems with sharp interfaces. Ray-splitting systems have recently attracted attention in the context of acoustic, quantum and electromagnetic wave chaos . It was shown that the mere presence of a ray-splitting boundary can drive an otherwise regular system into chaos . It was also shown that ray-splitting systems produce corrections to the Weyl formula  for the average density of states that can be computed analytically . The most conspicuous consequence of ray splitting is the existence of non-Newtonian periodic orbits in a ray-splitting system that contribute substantially to the fluctuating part of the level density . The existence of non-Newtonian orbits in a dielectric-loaded Bunimovich ray-splitting stadium was demonstrated experimentally . In addition it has been shown recently that exact trace formulae exist for a class of one-dimensional ray-splitting systems . These formulae were derived using results from quantum graph theory . Considering the two-point correlation function of the spectra of special quantum graphs, Schanz and Smilansky were able to derive novel combinatorial identities . This was possible by deriving the two-point correlation function in two independent ways, (i) directly using input from the quantum spectrum and (ii) using the exactness of the trace formula. Motivated by the methods of Kottos, Schanz and Smilansky  we show that novel combinatorial identities for binary Pólya necklaces  are obtained directly by comparing the spectral density of analytically solvable quantum graphs with their exact periodic orbit expansion.

The plan of this paper is as follows: in section II we present our model system, a one-dimensional ray-splitting system, whose spectrum can be obtained analytically. We present an exact periodic-orbit expansion of its level density expressed as a generalized Fourier sum over binary Pólya necklaces. In section III we outline our method for obtaining exact combinatorial identities for binary necklaces derived from analytically solvable cases of the one-dimensional ray-splitting system. In section IV we present two worked examples that yield two infinite sets of combinatorial identities. In section V we discuss our results and conclude the paper.

## II Spectrum

Denote by $`E_n`$ the spectrum of the one-dimensional scaled Schrödinger equation

$$\psi ^{\prime \prime }(x)+V_\lambda (E,x)\psi (x)=E\psi (x),$$
(1)
where

$$V_\lambda (E,x)=\{\begin{array}{cc}0,\hfill & \text{for }0<xa\text{,}\hfill \\ \lambda E,\hfill & \text{for }a<x<1\text{,}\hfill \\ \mathrm{},\hfill & \text{for }x(0,1)\hfill \end{array}$$
(2)
and $`\psi (x)0`$ for $`x(0,1)`$. Scaling potentials of the form (2) arise naturally in many ray-splitting systems, for instance in dielectric-loaded cavities . In this paper we concentrate on the case $`E>V_\lambda (E,x)`$ for all $`x`$. Then, without restriction of generality, the scaling constant $`\lambda `$ can be assumed to satisfy $`0\lambda <1`$. Define $`E=k^2`$, then the spectrum of (1) is determined by

$$\mathrm{sin}\left[k(\sigma _L+\sigma _R)\right]r\mathrm{sin}\left[k(\sigma _L\sigma _R)\right]=0,$$
(3)
where $`\sigma _L=a`$, $`\sigma _R=\eta (1a)`$, $`\eta =\sqrt{1\lambda }`$, and $`r=(1\eta )/(1+\eta )`$ is the reflection coefficient. For the derivations below it is useful to define the transmission coefficient $`t=\sqrt{1r^2}`$. Therefore, the reflection and transmission coefficients satisfy the relation

$$r^2+t^2=1.$$
(4)
Since in this paper we focus on the case $`0\lambda <1`$, both $`r`$ and $`t`$ are real and positive and range between 0 and 1. Possible quantum phases incurred at reflection or transmission events are treated explicitly and separately (see, e.g., (6) below). They are not included in $`r`$ or $`t`$.

We now discuss an alternative method of solving the quantum dynamics in the potential (2). This method is based on coding the periodic orbits with the help of symbol strings. We denote a bounce off $`x=0`$ by the letter $``$ and a bounce off $`x=1`$ by the letter $``$. Words formed with these two letters code for periodic orbits in (2). The word $``$, for instance, codes for the non-Newtonian orbit that bounces between $`x=0`$ and $`x=a`$ (above-barrier reflection orbit). The word $``$ codes for the (Newtonian) orbit that bounces between $`x=0`$ and $`x=1`$. Since a periodic orbit represented by the word $`w`$ cycles through the letters of $`w`$ without a well-defined beginning or end, two words $`w`$ and $`w^{}`$ are equivalent in our context, and code for the same periodic orbit, if they are of the same length (i.e. they consist of the same number of symbols) and their respective symbol sequences are identical up to cyclic permutations. Sequences of objects that are identical up to cyclic permutations are called (Pólya) necklaces . If the number of objects they consist of is two, they are called binary necklaces. Apparently therefore, the periodic orbits of (2) can be coded with the help of binary necklaces over the symbols $``$ and $``$. It is remarkable that for (2) every Newtonian or non-Newtonian periodic orbit can be mapped one-to-one onto a binary necklace. In other words, “pruning” is not necessary for the binary necklaces relevant for (2). Every binary necklace defines a possible periodic orbit of (2) and vice versa.

Given two letters, for instance $``$ and $``$, we can form $`2^{\mathrm{}}`$ words of length $`\mathrm{}`$. But, in general, many of these words will be cyclically equivalent, and correspond to the same necklace. So, how many necklaces of length $`\mathrm{}`$ are there? This question is answered by the following formula. There are exactly

$$N(\mathrm{})=\frac{1}{\mathrm{}}\underset{n|\mathrm{}}{}\varphi (n)\mathrm{\hspace{0.17em}2}^{\mathrm{}/n}$$
(5)
binary necklaces of length $`\mathrm{}`$, where the symbol “$`n|\mathrm{}`$” denotes “$`n`$ is a divisor of $`\mathrm{}`$”, and $`\varphi (n)`$ is Euler’s totient function defined as the number of positive integers smaller than $`n`$ and relatively prime to $`n`$ with $`\varphi (1)=1`$ as a useful convention. Thus the first four totients are given by $`\varphi (1)=1`$, $`\varphi (2)=1`$, $`\varphi (3)=2`$ and $`\varphi (4)=2`$. We illustrate (5) with two examples for $`\mathrm{}=1`$ and $`\mathrm{}=2`$. There are two necklaces for $`\mathrm{}=1`$, $``$ and $``$. Applying (5) to this problem, we verify $`N(1)=\varphi (1)\times 2=2`$. There are three necklaces of length 2, $``$, $``$ and $``$; again verified by (5), $`N(2)=[\varphi (1)\times 4+\varphi (2)\times 2]/2=3`$.

Given a binary necklace $`w`$, we define the following integer-valued functions on $`w`$: $`n_{}(w)`$ counts the number of $``$s in $`w`$, $`n_{}(w)`$ counts the number of $``$s, $`n(w)=n_{}(w)+n_{}(w)`$, $`\chi (w)`$ is the sum of $`n(w)`$ and the number of $``$-pairs in $`w`$, $`\alpha (w)`$ counts all occurrences of $``$-pairs or $``$-pairs, $`\beta (w)`$ counts all occurrences of $``$ or $``$ and $`\gamma (w)`$ is defined as $`\gamma (w)=2n_{}(w)+n_{}(w)`$. Note that the counting of $``$-pairs, $``$-pairs, $``$\- or $``$-combinations is to be understood cyclically, i.e., for example, $`\alpha ()=1`$ and $`\beta ()=2`$. Next we define the set $`W_p`$ of prime necklaces as the ones that cannot be written as a periodic concatenation of substrings. As shown recently , there exists an exact periodic orbit expansion for the spectral density of (1) in terms of prime binary necklaces,

$$\rho (k)=\overline{\rho }+\frac{1}{2\pi }\underset{wW_p}{}S_w\underset{\genfrac{}{}{0pt}{}{\nu =\mathrm{}}{\nu 0}}{\overset{\mathrm{}}{}}\left[(1)^{\chi (w)}r^{\alpha (w)}t^{\beta (w)}\right]^\nu e^{i\nu S_wk},k>0,$$
(6)
where

$$S_w=2[n_{}(w)\sigma _R+n_{}(w)\sigma _L]$$
(7)
is the action of the primitive periodic orbit coded by the prime binary necklace $`w`$ and $`\overline{\rho }=(\sigma _L+\sigma _R)/\pi `$ is the average level density.

## III Method

For special values of the parameters of the potential well (2) it is possible to solve (3) analytically, thus obtaining directly the density of states $`\rho (k)`$. Equating the explicit expression for $`\rho (k)`$ with the necklace expansion (6), one obtains combinatorial identities for binary necklaces. An illustrative example is the case $`a=1`$, for which the spectral density of (1) is given by

$$\rho (k)=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\delta \left(k\pi m\right).$$
(8)
In this case there exists only one primitive necklace, $``$, and the necklace expansion of (6) yields

$$\rho (k)=\frac{1}{\pi }\underset{\nu =\mathrm{}}{\overset{\mathrm{}}{}}e^{2i\nu k}.$$
(9)
Equating (8) and (9) yields the well-known Poisson formula.

A comment is in order here. For every finite $`a`$ the necklace expansion (6) involves an infinite sum over prime periodic necklaces. At $`a=1`$ this sum collapses to a single term. One may ask the question how this singular limit arises. The answer is the following. For every finite $`a`$ the actions of the right-hand lobes of the periodic orbits of (6) is finite. At $`a=1`$, these actions are zero. Consequently, all necklaces that represent different prime periodic orbits for $`a1`$ become repetitions of the Newtonian periodic orbit $``$ at $`a=1`$. This is the reason for the existence of only a single prime periodic orbit ($``$) at $`a=1`$.

Apart from trivial and well-known identities such as (8) and (9) above, (6) is a rich source of new and nontrivial combinatorial identities for binary necklaces. Specific examples are discussed in Section III. Here we outline the general method.

Equation (3) can be written as

$$\mathrm{sin}(\omega _1k)r\mathrm{sin}(\omega _2k)=0,$$
(10)
where $`\omega _1=\sigma _L+\sigma _R`$ and $`\omega _2=\sigma _L\sigma _R`$. A negative $`\omega _2`$ corresponds to a mirror reflection of $`V_\lambda (E,x)`$ with respect to $`x=1/2`$. Thus, because of $`\sigma _L`$,$`\sigma _R0`$, and without loss of generality, we may assume $`\omega _1\omega _20`$.

In case $`\omega _1`$ and $`\omega _2`$ are rationally related, i.e. $`\omega _1/\omega _2=p/q`$, $`pq𝐍`$ and $`p,q`$ relatively prime, (10) is reduced to the algebraic equation

$$\mathrm{sin}(p\omega k)r\mathrm{sin}(q\omega k)=0,$$
(11)
where $`\omega _1=p\omega `$ and $`\omega _2=q\omega `$. Using the formula

$$\mathrm{sin}(nx)=\mathrm{sin}(x)U_{n1}(\mathrm{cos}(x)),$$
(12)
where $`U_{n1}(x)`$ is the Chebyshev polynomial of the second kind, one obtains

$$\mathrm{sin}(\omega k)\left[U_{p1}(\mathrm{cos}\omega k)rU_{q1}(\mathrm{cos}\omega k)\right]=0.$$
(13)
It follows immediately from (13) that in the case of rationally related $`\omega _1`$ and $`\omega _2`$ there always exists a sequence of roots $`k_n^{(0)}=\pi n/\omega `$. The remaining roots are determined by

$$U_{p1}(x)rU_{q1}(x)=0,$$
(14)
where $`x=\mathrm{cos}(\omega k)`$. Since every root $`x_j`$ of (14) gives rise to a periodic sequence of eigenvalues, $`\mathrm{cos}(\omega k_n^{(j)})=x_j`$, $`j=1,2,\mathrm{},p1`$, together with the sequence $`k_n^{(0)}`$ the spectrum of (10) consists of $`p`$ (possibly degenerate) periodic sequences of roots. Whenever (14) can be solved analytically, the density of states

$$\rho (k)=\underset{j=0}{\overset{p1}{}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\delta \left(kk_n^{(j)}\right)$$
(15)
is known explicitly and together with (6) leads to a host of combinatorial identities for binary necklaces. Two examples are presented in the following section.

## IV Combinatorial identities

Example 1. For $`\sigma _L=\sigma _R`$ equation (3) becomes

$$\mathrm{sin}(2ka)=0,$$
(16)
with the solutions $`k_n=\pi n/(2a)`$. Note that there is no $`r`$-dependence in (16). The density of states is given by

$$\rho (k)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\delta \left(k\frac{\pi n}{2a}\right)=\frac{2a}{\pi }\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}e^{4imak}.$$
(17)
According to (7) and due to $`\sigma _L=\sigma _R`$, $`S_w`$ depends only on the binary length of $`w`$ and is given by $`S_w=2an(w)`$. Thus the sum (6) can be written as

$$\rho (k)=\frac{2a}{\pi }+\frac{a}{\pi }\underset{\genfrac{}{}{0pt}{}{m=\mathrm{}}{m0}}{\overset{\mathrm{}}{}}\underset{w\nu W_m}{}n(w)\left[(1)^{\chi (w)}r^{\alpha (w)}t^{\beta (w)}\right]^\nu e^{2imka},$$
(18)
where $`W_n`$ denotes the set of all length-$`n`$ binary necklaces w, $`w`$ is the shortest primitive code-piece in w and $`\nu `$ is the number of its repetitions in w. Comparing the series (17) and (18), we see that odd-length and even-length binary necklaces satisfy the sum rules

$$\underset{w\nu W_{2m1}}{}n(w)\left[(1)^{\chi (w)}r^{\alpha (w)}t^{\beta (w)}\right]^\nu =0,m=1,2,\mathrm{}$$
(19)
and

$$\frac{1}{2}\underset{w\nu W_{2m}}{}n(w)\left[(1)^{\chi (w)}r^{\alpha (w)}t^{\beta (w)}\right]^\nu =1,m=1,2,\mathrm{}.$$
(20)
At first glance it may seem surprising that (20) is a constant for all $`r`$. The solution lies in the relation (4). When properly ordered according to powers of $`r`$ and $`t`$, it turns out that (i) individual terms in (20) are of the form $`r^{2p}t^{2q}`$, where $`p+q=m`$ and (ii) the coefficients in front of the term $`r^{2p}t^{2q}`$ in (20) turn out to be binomial coefficients. Thus, for given $`m`$, the left-hand side of (20) reduces to $`(r^2+t^2)^m`$, which, according to (4) is equal to 1 for any choice of $`r`$. This explains why the seemingly variable left-hand side of (20) is nevertheless a constant. Thus we obtain from (20) the following infinite set of combinatorial identities for even-length binary necklaces

$$\frac{1}{2}\underset{w\nu W_{2m}}{}n(w)(1)^{\nu \chi (w)}\delta _{\nu \alpha (w)/2,s}=\left(\begin{array}{c}m\\ s\end{array}\right),s=0,\mathrm{},m,m=1,2,\mathrm{},$$
(21)
where $`\delta _{i,j}`$ is the Kronecker symbol.

In order to illustrate (21) let us first focus on the case $`m=1`$. According to (5) there are exactly three cyclically non-equivalent necklaces of binary length 2 given by $`\mathrm{w}_1=`$, $`\mathrm{w}_2=`$, $`\mathrm{w}_3=`$. The necklaces $`\mathrm{w}_1`$ and $`\mathrm{w}_3`$ are not primitive. The necklace $`\mathrm{w}_1`$ is a two-fold repetition of the primitive necklace $`w_1=`$. Thus $`\nu _1=2`$. An analogous consideration for $`\mathrm{w}_3`$ yields $`w_3=`$ and $`\nu _3=2`$. The necklace $`\mathrm{w}_2`$ is primitive. Therefore $`w_2=\mathrm{w}_2=`$ and $`\nu _2=1`$. The three necklaces $`\mathrm{w}_j`$, $`j=1,2,3`$, are listed in Table I. Also listed are their primitives $`w_j`$, the repetition indices $`\nu _j`$, and the values of the functions $`n(w_j)`$, $`\alpha (w_j)`$, $`\beta (w_j)`$, $`\gamma (w_j)`$ and $`\chi (w_j)`$.

We are now ready to check (21). For $`m=1`$ we have two choices for $`s`$: $`s=0`$ and $`s=1`$. For $`s=0`$ we have to scan the three words $`w_j`$, $`j=1,2,3`$, for $`\nu _j\alpha (w_j)/2=s=0`$. According to the entries in Table I, only $`\mathrm{w}_2`$ qualifies and the sum on the left-hand side of (21) reduces to the single term

$$\frac{1}{2}n(w_2)(1)^{\nu _2\chi (w_2)}=1=\left(\begin{array}{c}1\\ 0\end{array}\right).$$
(22)
This shows that (21) is indeed true for the simplest case $`m=1`$, $`s=0`$. For the case $`m=1`$, $`s=1`$ we have to check Table I for occurrences with $`\nu _j\alpha (w_j)/2=1`$. This is fulfilled for the necklaces $`\mathrm{w}_1`$ and $`\mathrm{w}_3`$. We obtain

$$\frac{1}{2}\left[n(w_1)(1)^{\nu _1\chi (w_1)}+n(w_3)(1)^{\nu _3\chi (w_3)}\right]=1=\left(\begin{array}{c}1\\ 1\end{array}\right).$$
(23)
This shows that (21) also works for $`m=1`$, $`s=1`$.

Testing (21) for $`m=2`$ involves finding all non-equivalent necklaces of binary length 4. According to (5) there are exactly six. All six necklaces are listed in Table II together with their properties. For $`m=2`$ we have three possibilities for $`s`$: $`s=0,1,2`$. For $`s=0`$ we have to check Table II for necklaces that fulfill $`\nu _j\alpha (w_j)/2=0`$. Only $`\mathrm{w}_4`$ qualifies. We obtain $`n(w_4)/2=1`$, which equals $`\left(\begin{array}{c}2\\ 0\end{array}\right)`$, the binomial coefficient on the right-hand side of (21). For $`s=1`$ we have to check for $`\nu _j\alpha (w_j)/2=1`$. We find three candidates: $`\mathrm{w}_2`$, $`\mathrm{w}_3`$ and $`\mathrm{w}_5`$. This time we have to be careful when summing the three terms on the left-hand side of (21), since $`\nu _3\chi (w_3)=5`$. Therefore the second term in (21) contributes with a minus sign. We obtain $`[n(w_2)n(w_3)+n(w_5)]/2=2`$, which equals $`\left(\begin{array}{c}2\\ 1\end{array}\right)`$, the corresponding binomial coefficient on the right-hand side of (21). Two necklaces, $`\mathrm{w}_1`$ and $`\mathrm{w}_6`$, contribute in the case $`s=2`$ and again satisfy (21).

Example 2. Suppose now that $`\sigma _L=2\sigma _R`$. In this case (3) becomes

$$\mathrm{sin}(ka/2)\left[4\mathrm{cos}^2(ka/2)r1\right]=0.$$
(24)
This equation has three sets of solutions

$$k_n^{(j)}=\frac{2j}{a}\mathrm{arccos}(\phi )+\frac{2\pi n}{a},j=1,0,1,$$
(25)
where $`\phi =\sqrt{1+r}/2`$. The density of states is

$`\rho (k)={\displaystyle \underset{j=1}{\overset{1}{}}}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}\delta \left(k+{\displaystyle \frac{2j}{a}}\mathrm{arccos}(\phi ){\displaystyle \frac{2\pi m}{a}}\right)=`$ (26)
$`{\displaystyle \frac{a}{2\pi }}{\displaystyle \underset{j=1}{\overset{1}{}}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e^{in[ak+2j\mathrm{arccos}(\phi )]}={\displaystyle \frac{a}{2\pi }}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e^{inak}\left[2T_{2n}(\phi )+1\right],`$ (27)
where $`T_n(x)\mathrm{cos}(n\mathrm{arccos}x)`$ are the Chebyshev polynomials of the first kind . Equating (27) order-by-order with the necklace expansion (6) we obtain the sum rules

$`{\displaystyle \underset{wW_p}{}}{\displaystyle \underset{\nu =1}{\overset{\mathrm{}}{}}}\gamma (w)\left[(1)^{\chi (w)}r^{\alpha (w)}(1r^2)^{\beta (w)/2}\right]^\nu \delta _{\nu \gamma (w),m}=`$ (28)
$`1+{\displaystyle \underset{j=0}{\overset{m}{}}}{\displaystyle \frac{2m(1)^j}{2mj}}\left(\begin{array}{c}2mj\\ j\end{array}\right)(1+r)^{mj},m=1,2,\mathrm{}.`$ (29)
We used formula 22:6:1 of ref.  for the Chebyshev polynomials in (27).

Ordering (29) according to powers of $`r`$, (29) can be reformulated as a combinatorial theorem on the set of binary necklaces, in which $``$ beads weigh twice as much as $``$ beads:

$$\underset{wW_p,C}{}\gamma (w)(1)^{[2m\chi (w)+s\gamma (w)m\alpha (w)]/[2\gamma (w)]}\left(\begin{array}{c}\frac{m\beta (w)}{2\gamma (w)}\\ \frac{s\gamma (w)m\alpha (w)}{2\gamma (w)}\end{array}\right)=$$

$$\delta _{s,0}+\underset{j=0}{\overset{ms}{}}\frac{2m(1)^j}{2mj}\left(\begin{array}{c}2mj\\ j\end{array}\right)\left(\begin{array}{c}mj\\ s\end{array}\right),s=0,1,\mathrm{}m,m=1,2,\mathrm{}.$$
(30)
The condition $`C`$ in the sum (30) is $`C=\gamma (w)|msm\alpha (w)/\gamma (w)`$ even. The sum on the left-hand side of (30) may be empty. In this case the sum is defined to be zero.

Let us check (30) with the help of a few examples. First we focus on the case $`m=1`$, $`s=0`$. In order to fulfill the first part of the condition $`C`$ in (30) we need $`\gamma =1`$. This, in turn, requires to find a necklace with $`n_{}=0`$ and $`n_{}=1`$. There is just one such necklace, namely $``$. But it does not fulfill the second part of $`C`$. Therefore, the sum on the left-hand side of (30) is empty, and the left-hand side is zero. The right-hand side adds up to $`1+12=0`$ and confirms (30) for this special case. For $`s=1`$ we find again that $``$ is the only choice for $`w`$. But this time the second part of $`C`$ is fulfilled and the left-hand side of (30) is

$$\gamma ()(1)^{[2\chi ()+\gamma ()\alpha ()]/[2\gamma ()]}\left(\begin{array}{c}\beta ()/[2\gamma ()]\\ [\gamma ()\alpha ()]/[2\gamma ()]\end{array}\right)=1.$$
(31)
We used $`\alpha ()=1`$, $`\beta ()=0`$, $`\gamma ()=1`$ and $`\chi ()=2`$. For $`m=1`$, $`s=1`$ the right-hand side of (30) consists of just one term, which turns out to be 1 as well. Thus we checked that (30) works for $`m=1`$. With the help of Tables I and II, other special cases may be checked as well.

## V Discussion and Conclusions

The work presented here is closely related to the theory of quantum graphs . While the quantum graphs considered by Kottos, Schanz and Smilansky  correspond to the case of zero potential on the bonds and delta potentials on the vertices, the step potentials considered in this paper correspond to constant potentials on the bonds and potential steps at the vertices. Thus, although the methods employed in this paper are essentially those used previously by Kottos, Schanz and Smilansky, we obtain a different class of combinatorial identities that apply to cyclic binary codes (Pólya necklaces). Another difference concerns the derivation of identities. While Kottos, Schanz and Smilansky use a route that involves two-point correlation functions, we show that novel combinatorial identities can be obtained directly from the periodic orbit expansions of explicitly solvable cases. These minor differences notwithstanding the central idea for generating entirely new classes of combinatorial identities is the same: Combinatorial identities can be obtained whenever a quantum system admits of (i) an explicit analytical solution and (ii) an exact periodic orbit expansion.

In addition to the two examples presented above, there exist many other cases in which (14) can be reduced to a low-order polynomial that can be solved by elementary means. Examples are the cases $`p=3`$, $`q=2`$ or $`p=5`$, $`q=3`$. Both cases can be treated in complete analogy to Example 2 above, and result in novel sum rules and combinatorial identities.

Recently we proved  that exact trace formulae exist for one-dimensional square wells with an arbitrary number of potential steps inside. Following the methods outlined above, our results can be generalized immediately to obtain novel combinatorial identities for necklaces with more than two types of beads.

## VI Acknowledgments

The authors acknowledge helpful comments and suggestions by Rick Jensen. Yu.D. and R.B. gratefully acknowledge financial support by NSF grants PHY-9900730 and PHY-9984075; Yu.D. by NSF grant PHY-9900746.

.

$$\begin{array}{ccccccccc}j& \mathrm{w}_j& w_j& \nu _j& n(w_j)& \alpha (w_j)& \beta (w_j)& \gamma (w_j)& \chi (w_j)\\ 1& & & 2& 1& 1& 0& 2& 1\\ 2& & & 1& 2& 0& 2& 3& 2\\ 3& & & 2& 1& 1& 0& 1& 2\end{array}$$

Table I: List of the three cyclically non-equivalent binary necklaces of length 2 together with their primitives ($`w`$) and repetition indices ($`\nu `$). Some properties of the primitives, such as their lengths ($`n`$), number of $``$ or $``$ pairs ($`\alpha `$), number of transmissions ($`\beta `$), their weighted lengths ($`\gamma `$) and their phases ($`\chi `$) are also listed.

$$\begin{array}{ccccccccc}j& \mathrm{w}_j& w_j& \nu _j& n(w_j)& \alpha (w_j)& \beta (w_j)& \gamma (w_j)& \chi (w_j)\\ 1& & & 4& 1& 1& 0& 2& 1\\ 2& & & 1& 4& 2& 2& 7& 4\\ 3& & & 1& 4& 2& 2& 6& 5\\ 4& & & 2& 2& 0& 2& 3& 2\\ 5& & & 1& 4& 2& 2& 5& 6\\ 6& & & 4& 1& 1& 0& 1& 2\end{array}$$

Table II: List of the six cyclically non-equivalent binary necklaces of length 4. The meaning of the columns is the same as in Table I.


text:
# Detection of H2 Emission from Mira B in UV Spectra from the Hubble Space TelescopeBased on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

## 1 INTRODUCTION

Mira (o Cet, HD 14386) is one of the most extensively observed variable stars, representing the prototype of the “Mira variable” class of pulsating stars. These stars are red giants on the asymptotic giant branch (AGB), with high mass loss rates and extended circumstellar envelopes. Mira’s pulsation properties are typical for long period Mira variables, with large luminosity changes of about 7 magnitudes during its 332 day pulsation cycle. Its distance is 128 pc according to Hipparcos (Perryman et al., 1997).

Mira has a hot companion star, Mira B (VZ Cet). The emission from Mira B is probably powered by a surrounding accretion disk fed by Mira’s massive wind. The accretor is generally assumed to be a white dwarf, although Jura & Helfand (1984) argue that the dearth of X-rays from the Mira system suggests that the companion is more likely a low mass main sequence star. The primary evidence for the accretion disk is very broad UV emission lines observed by the International Ultraviolet Explorer (IUE) (e.g., C IV $`\lambda `$1550, Si III\] $`\lambda `$1892, C III\] $`\lambda `$1909), which are believed to be from this rapidly rotating disk (Cassatella et al., 1985; Reimers & Cassatella, 1985).

Karovska et al. (1997) resolved the components of this binary at UV and optical wavelengths using HST Faint Object Camera (FOC) images, and also obtained the first spatially resolved spectra of each star individually. They measured a separation between the components in 1995 of only $`0.578^{\prime \prime }`$. The FOC images and the ultraviolet objective prism spectra provided for the first time detailed information on the continuum and line emission from each component of the system, and allowed a better measurement of the accretion luminosity of Mira B longward of 1500 Å. The resolution of the objective prism spectra was, however, insufficient for a detailed analysis of the line emission originating from the accretion disk.

In this paper, we present new high resolution spectroscopic observations of Mira B from the Hubble Space Telescope (HST) to further study the nature of accretion onto this intriguing star. In §2, we present the HST observations and compare them with previous UV spectra from IUE and FOC, and in §3 we discuss our results.

## 2 OBSERVATIONS AND RESULTS

### 2.1 Description of Observations

We observed Mira on 1999 August 2 using the Space Telescope Imaging Spectrograph (STIS) instrument on board HST. For a full description of STIS, see Kimble et al. (1998) and Woodgate et al. (1998). Unlike IUE, the STIS instrument has apertures small enough to observe the two components of the Mira AB system separately. Thus, there is no danger of any contamination by Mira A in our spectra. The Mira B observations consisted of two moderate resolution UV spectra taken through the $`0.2^{\prime \prime }\times 0.2^{\prime \prime }`$ aperture. The first utilized the E230M grating for a 600 s exposure of the 2303–3111 Å wavelength region. The second used the E140M grating for an 1850 s exposure of the 1140–1735 Å wavelength region. The spectra were processed using the standard HST/STIS pipeline software. To test the accuracy of our wavelength calibration, in the E140M spectrum we measured the centroid of the geocoronal H I Ly$`\alpha `$ emission peak to be at $`+27.8\pm 0.4`$ km s-1, in perfect agreement with its expected location for that time of year of $`+27.9`$ km s-1.

### 2.2 The H2 Lines

The top panel of Figure 1 shows the full E140M spectrum taken by HST/STIS. The spectrum is dominated by a very large number of narrow lines, which were not seen in previous IUE observations of Mira B. We identified them as being lines of molecular hydrogen. Furthermore, we believe the lines are all fluoresced by the H I Ly$`\alpha `$ line. For example, the 0-4 R(0) and 0-4 P(2) lines identified in Figure 1 are pumped by the 0-2 R(0) line at 1217.2046 Å, which lies within the very broad Ly$`\alpha `$ emission line from Mira B. Many other lines of the 0-x R(0) and 0-x P(2) sequences are also detected with flux ratios roughly consistent with the branching ratios (Abgrall et al., 1993). This provides strong support for the line identifications and for Ly$`\alpha `$ fluorescence being the source of the emission. Similarly, the 0-4 R(1) and 0-4 P(3) lines in Figure 1 are pumped by the 0-2 R(1) 1217.6437 Å line, and the 0-4 P(1) line is pumped by 0-2 P(1) at 1219.3677 Å.

Many Lyman band H2 lines were first identified in the solar spectrum. The solar H2 lines are fluoresced not only by Ly$`\alpha `$ but also by transition region lines (Jordan et al., 1977, 1978). Lines of H2 fluoresced by Ly$`\alpha `$ have by now been been detected in many other types of astrophysical objects, including red giant stars (McMurry et al., 1998; McMurry, Jordan, & Carpenter, 1999), T Tauri stars (Brown et al., 1981; Valenti, Johns-Krull, & Linsky, 2000), and Herbig-Haro objects (Schwartz, 1983; Curiel et al., 1995).

The H2 lines from Mira B all have widths of about 20 km s-1, well resolved by the $`7`$ km s-1 resolution of the E140M grating. They are centered at a heliocentric velocity of $`+57`$ km s-1. Various observations of Mira’s large, slowly expanding circumstellar envelope suggest a systemic radial velocity of about $`+56`$ km s-1 (Bowers & Knapp, 1988; Planesas et al., 1990; Josselin et al., 2000), consistent with our H2 velocity. A full analysis and list of the numerous H2 Lyman band lines in our FUV Mira B spectrum will be presented in a future paper.

### 2.3 Comparison with Previous Observations

The HST/STIS spectrum looks radically different from what we expected to see based on previous IUE observations (Cassatella et al., 1985; Reimers & Cassatella, 1985). Previous IUE observations of Mira B have instead shown a spectrum dominated by a few lines commonly observed in far-UV (FUV) spectra of stars and other astrophysical objects: C IV $`\lambda `$1550, O I $`\lambda `$1300, C II $`\lambda `$1336, etc. Figure 2 shows an example of one of the low resolution IUE spectra (SWP7029) of this spectral region, with line identifications based on previous analyses (Cassatella et al., 1985; Reimers & Cassatella, 1985).

In Figure 2, we rebinned and smoothed our STIS/E140M spectrum to match the resolution of the IUE spectrum. The deresolved STIS spectrum shows peaks at 1335 Å and 1400 Å, which Figure 1 shows to be predominantly from narrow H2 emission lines. In the IUE data, peaks at these locations were previously identified with C II and Si IV lines, respectively (with perhaps some O IV\] contribution to the 1400 Å peak). Were these features misidentified in the IUE spectra, and are they in fact mostly H2 emission in the IUE data as well as the STIS data? In order to address this question we computed another deresolved STIS spectrum for comparison with the IUE data, but only after removing all the narrow H2 lines. The result is shown as a solid line in Figure 2. The peaks at 1335 Å and 1400 Å decrease dramatically but do not disappear entirely, suggesting that there is some C II and Si IV emission in the STIS data. (The broad C II emission near 1337 Å can in fact be discerned in both panels of Fig. 1.) But what dominates these peaks in the IUE data, C II/Si IV or H2?

The most prominent H2 peak in the original STIS spectrum in Figure 2 (dotted line) that is not in danger of confusion with any other lines is the peak at 1275 Å, which completely disappears when the H2 lines are removed before deresolving the spectrum. There are 32 short-wavelength, low-resolution (SW-LO) spectra in the IUE archives. In none of these spectra is there a prominent peak at 1275 Å. Thus, we believe H2 emission is not an important contributor to any of the emission features observed by IUE, including the peaks at 1335 Å and 1400 Å. The STIS spectrum of Mira B is truly very different from anytime IUE observed it, and this is not just a consequence of the superior resolution and S/N properties of STIS compared with IUE.

Further evidence of this is provided by the continuum fluxes of our STIS spectra, which are well below those detected by IUE and previous HST observations. The fluxes of the IUE spectrum in Figure 2 were divided by 20 to match the continuum observed by STIS. The IUE data set clearly shows some variability (see Reimers & Cassatella, 1985), but only variations of about a factor of 2 or so. The STIS continuum fluxes are an order of magnitude lower than any of the IUE spectra. Figure 3 shows that this result extends to the near-UV (NUV) spectral region as well. The top panel of Figure 3 compares a rebinned, smoothed version of our STIS/E230M spectrum with a typical long-wavelength, low-resolution (LW-LO) IUE spectrum (LWP18360). Once again, the continuum fluxes are different by a factor of 20, and once again inspection of other LW-LO spectra in the IUE archive (31 total) shows this to be a general result. Comparison with the fluxes measured several years earlier with the HST also shows a significant drop in emission. The bottom panel of Figure 3 shows that the STIS continuum fluxes are a factor of 10 lower than those observed in 1995 by an HST Faint Object Camera (FOC) PRISM spectrum of Mira B presented by Karovska et al. (1997).

Emission line fluxes in our STIS spectra are also reduced relative to IUE, but to varying degrees. Most notable is the C IV $`\lambda `$1550 line, which dominates all the SW-LO IUE spectra, but is significantly weaker relative to the continuum in our STIS data. The C IV flux observed by STIS is roughly 100 times lower than observed by IUE (see Fig. 2). This behavior is different from that of cataclysmic variables, where C IV emission tends to become more prominent when the continuum decreases (Mauche & Raymond, 1987).

The optical emission from Mira B seems to also be low at the time of our STIS observations. In fact, we were not able to detect the companion in an image of Mira AB taken by HST using the narrow band \[O III\] filter of STIS centered at 5007 Å. We estimate a lower limit for the apparent magnitude of $`V>11`$. Mira B’s optical variations in the range $`V=912`$ have been suggested to have a possible period of 14 years (Joy, 1954; Yamashita & Maehara, 1977).

For a final comparison with IUE data, in Figure 4 we compare the Mg II k line profile observed by STIS (with a vacuum rest wavelength of 2796.352 Å) with a typical NUV, high resolution (LW-HI) observation from IUE (LWP29795). The two spectra are shown on a velocity scale centered on the rest frame of the star. The IUE fluxes are once again divided by 20 to roughly match the STIS fluxes. In the IUE profile, an opaque absorption feature is observed between 0 and $`400`$ km s-1 that is indicative of a high speed wind from Mira B. Similar absorption has been observed in optical H$`\beta `$ and H$`\gamma `$ lines (Yamashita & Maehara, 1977). Inspection of all 28 LW-HI IUE spectra available in the IUE archives reveals that all of the Mg II profiles have this absorption. The extent of the absorption does vary, sometimes extending to higher or lower velocities. However, the wind absorption is always opaque out to at least $`200`$ km s-1 in the IUE data. The STIS profile is different in that the absorption is only opaque to $`50`$ km s-1, and the absorption between 0 and $`50`$ km s-1 could conceivably be entirely interstellar. There is clearly wind absorption on the blue side of the STIS Mg II profile, but the wind opacity is much less than in any of the IUE spectra. This suggests that Mira B’s wind was much weaker at the time of the STIS observations than when IUE observed the star.

## 3 DISCUSSION

A fundamental question regarding the H2 lines that we have detected in the UV spectrum of Mira AB is, where are they coming from? The fluoresced H2 molecules may be located in the cool outer regions of the accretion disk surrounding Mira B, or perhaps the H2 is in Mira A’s slow, cool wind, which is fluoresced as it approaches the hot companion star and interacts with its faster wind (see below). Bowers & Knapp (1988) estimated that hydrogen in Mira’s wind should be roughly 70% molecular and 30% atomic. If the H2 is in Mira’s wind, one might expect to see a velocity shift of the H2 lines relative to the systemic velocity, which we do not observe. However, a shift is only expected if Mira B does not lie in the same plane of the sky as Mira A. Since Mira B’s orbit is not well known it is unclear what the geometry is. Furthermore, the velocity of Mira A’s wind is very slow, with a terminal expansion speed generally estimated to be in the $`27`$ km s-1 range (e.g., Ryde & Schöier, 2001). Thus, not much of a shift is to be expected and it remains very possible that the H2 lies within Mira’s wind.

The fluxes we observe in our HST/STIS observations of Mira B’s UV spectrum are dramatically lower than ever observed by IUE. Something may have happened to reduce the accretion rate onto the companion, thereby significantly lowering the accretion luminosity. The lower accretion rate leads to a lower mass loss rate from Mira B, based on the Mg II profiles (see Fig. 4).

It is not clear what could have caused the accretion rate to change. Instabilities in the accretion disk could be responsible. After all, more luminous accretion systems, such as dwarf novae and cataclysmic variables, show dramatic variability. It also might have been caused by inhomogeneities in the wind of Mira A that feeds the accretion disk, as there is evidence for substantial inhomogeneity in Mira A’s massive wind (e.g., Lopez et al., 1997).

A final fundamental question is why the lower accretion rate is also accompanied by more prominent emission from H2 lines. The lower accretion luminosity and weaker wind might allow more H2 molecules to exist close to Mira B, resulting in more fluorescence from Mira B’s broad Ly$`\alpha `$ emission. Future more detailed analysis and modeling of the very rich H2 fluorescence spectrum will shed some light on this issue.

We would like to thank J. Raymond for useful discussions and comments. Support for this work was provided by NASA through grant number GO-08298.01-99A from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. M. K. is a member of the Chandra Science Center, which is operated under contract NAS8-39073, and is partially supported by NASA.


text:
# Deconfinement transition and Luttinger to Fermi Liquid crossover in quasi one-dimensional systems

## Abstract

We investigate a system of one dimensional Hubbard chains of interacting fermions coupled by inter-chain hopping. Using a generalization of the Dynamical Mean Field Theory we study the deconfinement transition from a Mott insulator to a metal and the crossover between Luttinger and Fermi liquid phases. One-particle properties, local spin response and inter-chain optical conductivity are calculated. Possible applications to organic conductors are discussed.

The nature of the metallic phase of interacting electron systems depends strongly on dimensionality. In three dimensions, Fermi liquid (FL) theory applies, whereas in one dimension the quasi-particle concept breaks down, leading to a different kind of low-energy fixed point known as a Luttinger liquid (LL). For commensurate electron fillings, strong enough repulsive interactions destroy the metallic state altogether by opening a Mott gap. This phenomenon exists in all dimensions but the one-dimensional case is particularly favorable Voit (1995). In quasi one-dimensional (Q1D) systems, interchain hopping can induce a (deconfinement) transition from the Mott insulating (MI) state to a metallic state and crossovers between different metallic behaviors (Fig. 1). Besides their intrinsic theoretical interest, understanding these phenomena is directly relevant for a number of compounds such as the organic (super)-conductors (Bechgaard salts), which are three dimensional stacks of quarter- filled chainsJérome and Schulz (1982). Indeed, in these compounds some of the low temperature properties are well described by Fermi liquid theory, whereas optical Schwartz et al. (1998) and transport properties Moser et al. (1998) have shown that the high temperature phase is either a Luttinger liquid or a Mott insulator.

Describing these Q1D systems is not an easy task. The transverse hopping $`t_{}`$ is a relevant perturbation on the LL Bourbonnais et al. (1984); Wen (1990); Yakovenko (1992). Hence, perturbative renormalization group calculations yield an estimate of the crossover scale Bourbonnais et al. (1984); Bourbonnais and Caron (1991) but fail below that scale. In the MI state, electrons are confined on the chains by the Mott gap. A finite critical value of $`t_{}`$ is needed to induce an insulator to metal transition. Such a transition has been advocated Giamarchi (1997); Vescoli et al. (1998)to be at the root of the change from insulating to metallic behavior observed in Bechgaard salts when increasing pressure or when going from the TMTTF to the TMTSF family. Thus the effects of the inter-chain hopping that are the most important physically cannot be handled reliably by perturbative methods. Although some non-perturbative studies of deconfinement have been made for a finite number of chains Tsuchiizu et al. (2001) the case of an infinite system is still open. Some of the key questions yet to be answered are: (i) What is the crossover scale from the LL to MI, or from the LL to FL (Fig.1), and is there only one crossover scale for the different physical properties (transport, spin response, single-particle properties etc…) ? (ii) What is the nature of the low-temperature FL state, and in particular is the shape of the Fermi surface (FS) affected by interactions ? Do quasi-particle (QP) properties such as the quasi-particle residue $`Z(k)`$, vary significantly as the wavevector $`k`$ moves along the FS ?

In order to deal with these questions in a non-perturbative manner, a new method has been proposed Arrigoni (1999); Georges et al. (2000). This method, which generalizes the single site dynamical mean-field theory (DMFT) Georges et al. (1996), replaces the quasi one-dimensional system by a single effective chain from which electrons can hop to a self-consistent bath. We thus nickname it chain-DMFT (ch-DMFT ). It is exact in the limit of decoupled chains, while it reduces to usual DMFT in the opposite limit of zero-hopping along the chains, hence capturing the physics of a strongly correlated FL (or MI) Georges et al. (1996). Solving the mean field equations is a very challenging task however. In this paper, the first quantitative solution of the ch-DMFT equations is presented. We show that the above issues can be valuably addressed, and some of the above questions answered. In view of the difficulty of the problem only the Hubbard model is considered in this paper, but more realistic extensions are in sight, as discussed below.

Let us consider a system of coupled chains described by a sum of Hamiltonians of the isolated chains plus a hopping term between neighboring chains. The ch-DMFT approximates this system by the effective action:

$`S_{\text{eff}}=`$ $`{\displaystyle _0^\beta 𝑑\tau 𝑑\tau ^{}\underset{ij,\sigma }{}c_{i\sigma }^+(\tau )𝒢_0^1(ij,\tau \tau ^{})c_{j\sigma }(\tau ^{})}`$ (1)
$`+{\displaystyle _0^\beta }𝑑\tau H_{1D}^{int}[\{c_{i\sigma },c_{i\sigma }^+\}]`$
where $`H_{1D}^{int}`$ is the interacting part of the on-chain Hamiltonian. $`𝒢_0`$ is an effective propagator determined from a self-consistency condition: the Green’s function $`G(ij,\tau \tau ^{})c(i,\tau )c^+(j,\tau ^{})_{\text{eff}}`$ calculated from $`S_{\text{eff}}`$ should coincide with the on-chain Green’s function of the original problem, with the same self-energy $`\mathrm{\Sigma }=𝒢_0^1G^1`$. This reads:

$$G(k,i\omega _n)=𝑑ϵ_{}\frac{D(ϵ_{})}{i\omega _n+\mu ϵ_k\mathrm{\Sigma }(k,i\omega _n)ϵ_{}}$$
(2)
Let $`ϵ_{}(k_{})`$ denote the Fourier transform of the inter-chain hopping $`t_{}^{mm^{}}`$, $`D(ϵ_{})=_k\delta [ϵ_{}ϵ_{}(k_{})]`$ the corresponding density of states, $`k`$ the momentum in the chain direction and $`\omega _n`$ the Matsubara frequencies corresponding to the inverse temperature $`\beta `$. The ch-DMFT equations (1,2) fully determine the self-energy and Green’s function of the coupled chains $`G(k,k_{},i\omega _n)^1=i\omega _n+\mu ϵ_kϵ_{}(k_{})\mathrm{\Sigma }(k,i\omega _n)`$. The ch-DMFT approach can be rigorously justified in the formal limit where the lattice connectivity in the transverse direction is taken to infinity. In systems with finite transverse connectivity it can be viewed as an approximation neglecting the dependence of the self-energy on transverse momentum, keeping both frequency and in-chain momentum dependence: $`\mathrm{\Sigma }=\mathrm{\Sigma }(k,i\omega _n)`$. In the following, we specialize to the Hubbard model, i.e for each chain:

$$H_{1D}=t\underset{i}{}\left[c_i^+c_{i+1}+c_{i+1}^+c_i\right]+U\underset{i}{}n_in_i$$
(3)
and we consider a two-dimensional array of chains forming a square lattice with nearest neighbor hopping, so that the transverse dispersion is $`ϵ_{}(k_{})=2t_{}\mathrm{cos}k_{}`$.

We solve the effective chain problem for chains of 16 or 32 sites with periodic boundary conditions using the QMC Hirsch-Fye algorithm Hirsch and Fye (1986). Using 32 time-slices in imaginary time allows to reliably access temperatures down to $`T/W1/50=0.02`$ with $`U/W1`$ where $`W=4t`$ is the bandwidth of the one-dimensional chain. Typically, about 5000 QMC sweeps and 10 to 15 ch-DMFT iterations are sufficient to reach convergence. The main quantities that we calculate and analyze are: a) On-chain (i.e summed over $`k_{}`$) single particle Green’s functions: $`G(k,i\omega _n)`$, and self-energies $`\mathrm{\Sigma }(k,i\omega _n)`$ (obtained from $`𝒢_0^1G^1`$ 111To get a reliable estimate for the self-energy we used an analytical Fourier transform of a spline interpolated Green’s function with exact boundary conditions in the spirit of Joo et al., condmat-0009367.). This allows to identify the location of the FS $`k_{}^F(k)`$ by solving $`ϵ_{}(k_{}^F)=\mu \mathrm{Re}\mathrm{\Sigma }(k,i\omega _{n=1})ϵ_k`$ and the QP residue $`Z=Z(k_{}^F)`$ by fitting the slope ($`=1Z^1`$) of $`\mathrm{\Sigma }`$ vs. $`i\omega _n`$ at a specific FS point. b) On-chain spin and charge response functions and in particular the local spin response function: $`\chi _s(\tau )=S^z(j,0)S^z(j,\tau )=_{k,k}\chi _s(k,k_{},\tau )`$. Indeed, in a LL (or with $`K_\rho =1`$ in a FL) we have $`\chi _s(\tau )=\chi _s(\beta /2)\left(\mathrm{sin}\pi \tau /\beta \right)^{(1+K_\rho )}`$ in the asymptotic regime where $`\beta ,\tau `$ are larger than some coherence scale. By fitting our results to this form, we can define an effective $`K_\rho `$. Furthermore, the temperature dependence of $`\chi _s(\beta /2)`$ is closely related to the NMR relaxation rate $`1/T_1lim_{\omega 0}\chi _s^{\prime \prime }(\omega )/\omega `$. In a FL liquid at low-enough temperature Randeria et al. (1992): $`1/T_1=\beta \chi _s(\beta /2)/2\pi ^2`$. In a LL $`1/T_1`$ and $`\beta \chi _s(\beta /2)`$ have the same T-dependence (albeit with a proportionality factor depending on $`K_\rho `$). c) Inter-chain optical conductivity. Within ch-DMFT vertex corrections drop out Georges et al. (2000) and the conductivity is given by

$`\sigma _{}(i\omega )`$ $`{\displaystyle \frac{t_{}^2}{\omega }}{\displaystyle \frac{dk_{}}{2\pi }\mathrm{sin}^2k_{}\frac{dk}{2\pi }\frac{1}{\beta }\underset{m}{}G(k,k_{},i\omega _m)}`$ (4)
$`\times G(k,k_{},i\omega _m+i\omega )`$
Note that we have taken into account the $`k_{}`$\- dependence of the current vertex in (4). In practice, we perform analytically the $`k_{}`$-integration in (4) and numerically the k-integration and Matsubara summation. We then use the Maximum Entropy (MaxEnt) algorithm Jarrell and Gubernatis (1996) to continue $`\sigma _{}(i\omega )`$ to the real axis and obtain the inter-chain optical conductivity $`\mathrm{Re}\sigma _{}(\omega ,T)`$.

We first discuss our results away from half-filling, for $`U/W=1`$ and $`\mu =0.2`$ (corresponding to a total density $`n0.8`$). The numerical method was thoroughly tested for decoupled chains ($`t_{}=0`$), for which our fitting procedure of $`\chi _s`$ yields $`K_\rho 0.7`$, in agreement with the exact result Schulz (1990). We then study how inter-chain coherence develops as temperature is lowered for coupled chains with $`t_{}/W=0.14`$. Fig. 2 displays the effective $`K_\rho `$ as a function of temperature: a clear crossover from a LL (with $`K_\rho 0.7`$) to a FL ($`K_\rho =1`$) is seen as temperature is lowered below $`T^{}/W1/44`$.

According to Bourbonnais et al. (1984); Bourbonnais and Caron (1991), the crossover scale is given by $`T^{}=\frac{t_{}}{\pi }C(t_{}/t)^{\theta /(1\theta )}`$, with $`\theta =(K_\rho +1/K_\rho )/41/2`$ ($`0.03`$ here). The interactions reduce its value compared to the naive one $`t_{}/\pi `$. We cannot meaningfully test this renormalization effect because of the small values of $`\theta `$ in the Hubbard model. Our results are consistent with $`T^{}Ct_{}/\pi `$ with a proportionality constant $`C0.5`$. In the low temperature FL regime, we find that the location of the FS is essentially unaffected by interactions. The quasi-particle residue $`Z`$ is not uniform along the FS. Its dependence on $`k_{}^F`$ is displayed in Table 1. $`Z`$ is larger for $`k_{}`$ close to $`\pi `$, while it is somewhat smaller for small $`k_{}`$. This is not a very big effect however, and it would be inappropriate to speak of “hot spots” Zheleznyak and Yakovenko (1995) in that region (even though it may be indicative that such a phenomenon appears for stronger interactions). Note that the single-chain approximation $`\mathrm{\Sigma }=\mathrm{\Sigma }_{1D}`$ used by many authors would lead to a much stronger variation of $`Z`$ along the FS. Our results, in agreement with those of Arrigoni (1999), show that $`Z`$ is more uniform than this approximation would suggest.

We now turn to the commensurate (half-filled) case, and present results for $`U=0.65W`$. By fitting the self-energy to an analytical form, we first checked that our numerical results are consistent with a Mott insulating state with a gap $`\mathrm{\Delta }_{1D}/W=0.1`$, in agreement with exact results. Turning on $`t_{}`$, we estimate the $`t_{}`$-dependence of the gap by performing simulations at one of the lowest accessible temperatures $`T/W=1/40`$. The gap vanishes for $`t_{}/W0.07`$. Hence our method captures the insulator-to-metal transition induced by transverse hopping.

The effective $`K_\rho `$, shown in Table 2 is also a clear indicator of the MI ($`K_\rho =0`$) to FL ($`K_\rho =1`$) transition, in agreement with the qualitative expectations of Fig. 1. At this low temperature, the intermediate LL regime is too narrow to be visible. The location of the transition is in reasonable agreement with the naive criterion $`\mathrm{\Delta }_{1D}t_{}^{}`$ Giamarchi (1997); Vescoli et al. (1998); Tsuchiizu et al. (2001). Fig. 3 displays the inter-chain optical conductivity for several values of $`t_{}`$. In the MI phase, $`\sigma _{}(\omega )`$ shows a gap, followed by an onset of absorption starting at approximately the gap and extending up to a scale of order $`U`$, where a broad second peak is apparent. A low-frequency Drude peak develops as the insulator to metal transition is crossed. Close to the transition, the weight in the Drude peak is small, while the Hubbard band feature is still clearly visible and carries a significant part of the spectral weight.

The results for the NMR relaxation rate $`1/T_1`$ are displayed in Fig. 4 for a value of $`t_{}/W=0.11`$.

At low temperatures the Korringa law $`1/(T_1T)=const`$ is recovered, while at higher temperatures $`1/T_1T`$ is temperature dependent. The observed T-dependence is consistent with $`1/(T_1T)T^{K_\rho 1}`$ Bourbonnais et al. (1984) with a T-dependent exponent $`K_\rho 1`$, indicative of a LL evolving gradually into a FL (cf. Fig. 1). As for the doped case, the FS is nearly indistinguishable from the non-interacting case (see Fig. 5). The QP residue $`Z_k_{}`$ depends only very weakly on the location on the Fermi surface as shown in Table 3.

The data may suggest very shallow minima, which in contrast to the doped case, are found at $`k_{}\pm \pi /2`$, corresponding to a vanishing inter-chain kinetic energy. We would like to contrast our findings with the conclusions drawn from the “single-chain” (RPA) approximation Wen (1990); Essler and Tsvelik (2001): $`\mathrm{\Sigma }=\mathrm{\Sigma }_{1D}`$. As recently pointed out Essler and Tsvelik (2001), this approximation predicts that the FS close to the metal-insulator transition should consist of disconnected electron and hole “pockets”, as depicted schematically on Fig. 5. This is because the self-energy of a MI diverges at low-frequency for $`k=\pi /2`$, so that no FS point can correspond to this value of $`k`$. In contrast, our data find a conventional connected FS down to the transition point (up to our numerical accuracy). This shows that it is crucial, as done in ch-DMFT , to take into account the feedback effects of the interchain hopping in the self-energy.

In conclusion, we have shown in this paper that the ch-DMFT approach is a tool of choice for the study of the crossovers and insulator to metal transitions induced by transverse hopping in quasi one-dimensional systems. Several features found in our study are qualitatively reminiscent of experimental observations on organic conductors, such as: the deconfinement transition itself, the crossover from a LL at high-T to a FL at low-T (as revealed e.g in NMR) and most significantly the coexistence of a Drude feature with small spectral weight and of a Hubbard-band feature in optical conductivity. Extensions of the model including on-site and nearest neighbor interactions would allow to study the deconfinement transition at quarter-filling and thus to make a more realistic Giamarchi (1997); Schwartz et al. (1998) comparison with experimental results. Though this is computationally much more demanding, we hope to use ch-DMFT in this context in future work.

###### Acknowledgements.

Acknowledgements: We acknowledge useful discussions with A. Tsvelik, F. Essler, D. Jerome and C. Bourbonnais. This research has been supported by a Marie Curie Fellowship of the European Community Programme “Improving Human Potential” under contract number MCFI-2000804 and a grant of supercomputing time at NIC Jülich.