text:
# X-Ray Lines from Hot Flows around White Dwarfs. Application to SS Cygni.

## 1 Introduction

Dwarf novae (DN) are members of the class of Cataclysmic Variables (CVs), which are binary stars composed of a main-sequence donor which transfers mass via Roche-lobe overflow onto a White Dwarf (WD). DN experience quasi-regular, luminous outbursts during which accretion onto the WD proceeds at a high rate. Most of the time, however, DN are in quiescence, a phase during which the accretion rate onto the WD is much reduced (see Warner 1995 for a review).

Quiescent DN are well known sources of hard X-ray emission (Córdova & Mason 1983; Patterson & Raymond 1985a). This emission, typically at a level of $`10^{3032}`$ erg s-1, is consistent with Bremsstrahlung emission from a hot gas with temperatures $`220`$ keV (Patterson & Raymond 1985a; Eracleous, Halpern & Patterson 1991; Belloni et al. 1991; Yoshida, Inoue & Osaki 1992; Mukai & Shiokawa 1993). It has been commonly attributed to the boundary layer (BL) of the accreting WD, because at low enough accretion rates ($`<\mathrm{\hspace{0.33em}10}^{16}`$ g s-1), the BL is unable to cool before it expands out of the plane of the accretion disk, so that it is hot, optically-thin and therefore a significant source of hard X-rays (Pringle & Savonije 1979; Tylenda 1981; King & Shaviv 1984; Patterson & Raymond 1985a).

Detailed calculations by Narayan & Popham (1993) show that the optically-thin BLs of accreting WDs are also radially extended (on the order of the WD radius in their models) and that they advect part of the internally dissipated energy inwards to the WD surface as a consequence of their inability to cool. The similarities between the optically-thin BL solutions of Narayan & Popham (1993) and ADAF solutions (Ichimaru 1977; Rees et al. 1982; Narayan & Yi 1994; 1995; Abramowicz et al. 1995; see Narayan, Mahadevan & Quataert 1998 for a review) prompted Menou (2000) to propose that ADAFs could actually be present around WDs accreting at low rates, such as in some quiescent DN.

ADAF models applied to accreting WDs are able to reproduce the typical X-ray luminosities of quiescent DN, as well as the range of plasma temperatures inferred from early X-ray observations. In addition, advection in the ADAF results in a large amount of energy deposited onto the central WD. If this energy is thermalized in the WD atmosphere before being radiated away, a strong EUV emission component originates from the heated atmosphere which can explain the luminous He II $`\lambda 4686`$ lines often observed in quiescent dwarf novae as being the result of disk irradiation (Menou 2000). In fact, the presence of such a component had been suggested by previous studies which emphasized the large supply of EUV photons (difficult to observe directly) required to explain the strength of the lines by photoionization of the disk material (Patterson & Raymond 1985b; Vrtilek et al. 1994). Simple BL models (in which one assumes the kinetic energy of the accreting gas just above the WD surface is radiated as isothermal, optically-thin Bremsstrahlung emission; e.g., Tylenda 1981; Patterson & Raymond 1985a) do not offer any satisfying explanation to the strength of these lines (see also, e.g., Ferland et al. (1982) for a discussion of additional difficulties faced by simple models).

The success of ADAF models applied to quiescent DN in providing an explanation for both the hard X-ray emission and the strong He II lines is attractive, but not compelling. In particular, X-ray eclipse observations indicate that a radially extended ADAF is not present in all quiescent DN (Wood et al. 1995; Mukai et al. 1997; Pratt et al. 1999; Menou 2000). In this paper, we emphasize that the presence of an ADAF can be tested in non-eclipsing systems by detecting (or not) thermal emission lines from the hot plasma, which are robust predictions of the theoretical models (Narayan & Raymond 1999). Line diagnostics can put severe constraints on the structure of the hot flow (Perna, Raymond & Narayan 2000); this is especially true for quiescent DN, which are systems with a value of the WD mass and the inclination to the line of sight often relatively well constrained.

Other variations on the theme of a thin Shakura-Sunyaev (1973) disk with a hot, optically thin boundary layer include the model of King & Shaviv (1984), who emphasized the role of thermal conduction and the qualitative change expected at higher accretion rates. The optically thick boundary layers of high accretion rate systems produce less luminosity than expected from the outer disk luminosity (e.g., Mauche, Raymond & Mattei 1995; Popham & Narayan 1995) through some combination of advection and loss in a wind.

In §2, we recall some important properties of the DN SS Cyg to be used as constraints for the models. In §3, we describe the ADAF models and the numerical code used for the spectral predictions. In §4, we show predictions for the system SS Cyg and we emphasize the important constraints that can be put on the structure of the hot flow from detailed X-ray line diagnostics. Finally, we discuss some consequences and limitations of this work in §5, before summarizing our main results in §6.

## 2 SS Cygni and its quiescent X-ray emission

Observations show that some quiescent DN are more likely to harbor an ADAF than others. In particular, of the five DN observed in quiescence by the Hopkins Ultraviolet Telescope (Long 1996; spectral range $`8301860\AA `$), three (WX Hyi, Yz Cnc and SS Cyg) show strong, blue continua with no evidence for the presence of a normal WD atmosphere (broad instead of narrow emission lines and no evidence for the expected broad Ly$`\alpha `$ absorption feature). These three sources are ideal candidates for the presence of an ADAF because a blue UV continuum is expected if the WD atmosphere is heated up to several $`10,000`$ K by a substantial amount of energy advection in the surrounding hot flow (Menou 2000). SS Cyg is one of the best studied DN (see, e.g., Warner 1995), a bright source of X-ray emission in quiescence and is an approved target for observation with the Chandra X-ray Observatory. This motivates us to construct specific models for this system in what follows. The only other X-ray spectral models of SS Cyg in quiescence that we know of are the one-temperature models of Done & Osborne (1997), which are based on the BL solutions of Narayan & Popham (1993).

The mass of the WD in SS Cyg is $`M_{\mathrm{WD}}1.2M_{}`$ and the inclination of the system to the line of sight is $`i40^o`$ (Ritter & Kolb 1998). This inclination guarantees that the broadening of the X-ray lines from an ADAF is dominated by the bulk motion of the flow rather than thermal and turbulent motions (see below). The X-ray emission of SS Cyg during quiescence is variable. It has been characterized by an X-ray luminosity in the range $`10^{31}10^{32}`$ erg s-1 (in various X-ray bands) and gas temperatures in the range $`519`$ keV (see, e.g., Yoshida et al. 1992; Mukai & Shiokawa 1993).

## 3 Modeling techniques

There are a number of important assumptions made in the models presented below. For clarity, we briefly recall them here, while a more detailed discussion of these assumptions and their validity can be found in Menou (2000).

The hot flow structure used in the present study for the spectral calculations is nearly identical to the self-similar ADAF solutions described by Narayan & Yi (1994). This flow must make a transition from a large (but still subsonic) radial speed to a zero radial speed at the stellar surface. It is assumed that this transition happens in a narrow region (boundary layer, or BL) in the vicinity of the WD surface. The structure of this BL is not described in any detail here. The amount of mechanical energy released in this region (which depends on the stellar rotation rate) is comparable to the amount of energy radiated by the ADAF (Menou 2000). In all the models presented, we assume that both this mechanical energy and the much larger amount of energy advected by the ADAF (see below) are released at the WD surface as extreme-UV radiation,111The large amount of energy advection results in significant heating of the WD and could have important consequences on the structure of its envelope (e.g. Shaviv & Starrfield 1987; Pringle 1988; Regev & Shara 1989) so that these additional components are not considered in the X-ray emission models presented (broadband models show that the EUV radiation goes to infinity without interacting significantly with the flow, though prominent HeII emission lines may be a signature of this component; Menou 2000). We expect a roughly blackbody spectrum with the luminosity of the advected energy from the WD surface, but detailed models of the fraction of the WD surface covered and of the structure of the heated atmosphere are needed to make more specific predictions.

### 3.1 ADAF models

The radial profiles of physical quantities in the hot flow, such as density, radial and azimuthal velocities and temperature are obtained from ADAF models and then used as an input for the detailed spectral calculations described in §3.2. We use two specific ADAF models presented by Menou (2000), which have parameters appropriate for SS Cyg: $`M_{\mathrm{WD}}=1.2M_{}`$ (WD mass), $`R_{\mathrm{WD}}=5\times 10^8`$ cm (WD radius), $`\alpha _{\mathrm{ADAF}}=0.2`$ (viscosity parameter), $`\beta =6`$ (ratio of gas to magnetic pressure), $`\gamma =1.636`$ (adiabatic index of the fluid, which includes contributions from the particles, the turbulence and the magnetic field; see Quataert & Narayan 1999), $`\delta =10^2`$ (fraction of direct electron viscous heating; unimportant for the present models), $`p=0`$ (no wind) and $`\dot{M}=7\times 10^3\dot{M}_{\mathrm{Edd}}`$ (accretion rate), where $`\dot{M}_{\mathrm{Edd}}=1.39\times 10^{18}(M_{\mathrm{WD}}/M_{})`$ g s-1 is the Eddington accretion rate for a $`10\%`$ radiative efficiency. In the two models, this value of $`\dot{m}`$ corresponds to a physical accretion rate $`\dot{M}10^{16}`$ g s-1 $`1.5\times 10^{10}M_{}`$ yr-1, a $`0.510`$ keV X-ray luminosity of $`2\times 10^{31}`$ erg s-1 and $`85`$% of the total gravitational energy liberated being advected by the flow.

The difference between the two models comes from the radial extent of the ADAF. In model I, the ADAF extends from the WD surface at $`R_{\mathrm{WD}}`$ up to $`R_{\mathrm{max}}=10^4R_s`$, where $`R_s=3\times 10^5(M_{\mathrm{WD}}/M_{})`$ cm is the Schwarzschild radius. In model II, the ADAF extends further out, up to $`R_{\mathrm{max}}=10^5R_s`$. In addition, the accretion rate in the ADAF decreases as $`R_{\mathrm{trans}}/R`$ beyond $`R_{\mathrm{trans}}=10^4R_s`$ in model II, to represent a corona accreting gas that is gradually being evaporated from an underlying thin disk (see Menou 2000 and Esin, McClintock & Narayan 1997 for details). In the following, we neglect the role, even indirect, that the disk could have on the X-ray spectrum emitted by the system (see §5 for a discussion of a possible reflection component). Except for the spectral calculations, all the other model characteristics are similar to those described in Menou (2000).

Panels (a) and (b) of Figure 1 show the radial profiles of various quantities in models I and II, respectively. At radii $`10^3R_s`$, the ions and the electrons in the ADAF have essentially the same temperature (because of an efficient coupling by Coulomb collisions), from a few $`10^8`$ K close to the WD to $`10^7`$ K further away. This range of temperatures is ideal for the efficient production of thermal X-ray emission lines (Narayan & Raymond 1999). The radial profiles closely follow the self-similar ADAF solutions of Narayan & Yi (1994). A change in the slope of the density profile is clearly visible in Figure 1b; it corresponds to the radius $`R_{\mathrm{trans}}`$ beyond which the mass accretion rate in the ADAF starts decreasing with $`R`$. Figure 1 also shows that the azimuthal speed of the gas in the ADAF (or equivalently its angular rotation speed) is substantially sub-Keplerian ($`\mathrm{\Omega }0.2\mathrm{\Omega }_K`$ in the present models).

### 3.2 Spectral models

We compute the X-ray spectra using an extended version of the Raymond & Smith (1977) code. The code computes the Bremsstrahlung, recombination and two-photon continua, as well as the emission in spectral lines. The effects of collisional excitations (see, e.g. Pradhan, Norcross & Hummer 1981 and Pradhan 1985), recombinations to excited levels of H-like and He-like ions (see, e.g., Mewe, Schrijver, & Sylwester 1980) and dielectronic recombination satellite lines (see, e.g. Dubau et al. 1981 and Bely-Daubau et al. 1982) are all included in the calculation. The ionization state of the elements is computed self-consistently within the code itself. Given the typical range of densities and infall velocities in the ADAF (see Fig. 1), the ionization time is much smaller than the infall time and the ionization state of the elements is very close to equilibrium. Note that in computing the expected luminosities, we do not take into account the partial obscuration by the central WD of the inner shells of the accretion flow. This should not strongly affect our predictions because the outer regions of the flow contribute significantly to the total luminosity in the radially extended ADAF models presented here. The high temperature regions close to the WD surface would be most strongly affected and the fluxes from these regions (which may be so hot that they do not contribute much to the line emission) would be reduced by a factor two at most.

Let $`E_m(r,\theta )`$ be the emissivity in a given line, and let $`v_{\mathrm{los}}`$ be the local component of the velocity along the line of sight to the observer. The emission profile of a line as measured by an observer is then given by

$$\mathrm{\Phi }(v)=_0^{2\pi }𝑑\varphi _0^\pi 𝑑\theta \mathrm{sin}\theta _{r_{\mathrm{in}}}^{r_{\mathrm{out}}}𝑑rr^2E_m(r,\theta )\mathrm{\Phi }^{}(v;r,\theta ,\varphi )\mathrm{ergs}\mathrm{s}^1(\mathrm{km}\mathrm{s}^1)^1,$$
(1)
where

$$\mathrm{\Phi }^{}(v;r,\theta ,\varphi )=\frac{1}{\sqrt{\pi }\mathrm{\Delta }}\mathrm{exp}\left\{\frac{[vv_{\mathrm{los}}(r,\theta ,\varphi )]^2}{\mathrm{\Delta }^2}\right\},$$
(2)
and $`\mathrm{\Delta }=\sqrt{2k_\mathrm{B}T/m_a+v_{\mathrm{turb}}^2}`$ for an atom of mass $`m_a`$. For the turbulent velocity, we adopt $`v_{\mathrm{turb}}\alpha _{\mathrm{ADAF}}c_s`$, where $`c_s=(\gamma \mu k_BT/m_p)^{1/2}`$ is the fluid adiabatic sound speed and $`\mu =0.5`$ its mean molecular weight. In the models presented below, we find that line broadening is largely dominated by the flow bulk motion, with only a small contribution from the turbulent motions (the contribution from the thermal motions is negligible for the heavy element X-ray lines discussed here; this also implies that turbulence is supersonic for the heaviest species).

## 4 Results

### 4.1 Continuum and Line Emission

Whereas the emission of thermal X-ray lines is a robust prediction of models of hot accretion flows, the detailed spectrum is not. It depends on several parameters of the model, such as the radial extent and the run of density in the hot flow. As we show below, this property can be used, in turn, to constrain the structure of the hot flow from detailed X-ray spectroscopic diagnostics.

The X-ray spectra corresponding to models I and II are shown in panels (a) and (b) of Figure 2 . The spectra are characterized by emission lines of Fe, O, Si and other heavy elements, superposed on the Bremsstrahlung continuum. Note that a solar abundance has been assumed for the elements in the gas. Deviations from solar abundances could result in the lines of specific elements being weaker or stronger than in the models presented here. The shape of the continuum and line emission contains information on the structure of the hot flow. As a sum of the emission of shells of various densities $`\rho `$ and temperatures $`T`$ (see Fig. 1), a careful fit to an observational spectrum should provide valuable constraints on $`\rho `$ and $`T`$ in the hot flow. For instance, it is clear from Fig. 2 that there is additional low-energy Bremsstrahlung continuum emission in model II relative to model I, because of the additional emission from the outer and cooler regions of the flow.

The relative importance of the various lines in the spectrum is also indicative of the density, temperature and radial extent of the hot flow. For the ADAF models, lines formed at temperatures $`10^8`$ K are absent from the emission spectrum if the radial extent of the flow is $`10^4R_s`$ (Fig. 2, panel a), while they become clearly apparent if the outer radius of the flow corresponds to a temperature $`10^7`$ K (model shown in Fig. 2, panel b). Similarly, in a model where the accretion rate in the ADAF is not decreasing with radius beyond $`R_{\mathrm{trans}}=10^4R_s`$ (like in model II), but is constant with radius, lines preferentially emitted at temperatures $`10^8`$ K would be stronger because of the higher densities in the outermost regions of the flow. This illustrates that the relative strengths of lines formed at various radii can be used to constrain the radial extent of the flow and its density profile222 Note that this can be done even for unknown elemental abundances, by considering lines from different ions of the same element. (see Perna, Raymond & Narayan 2000 for details).

This can be complemented by measurements of line profiles, which strongly constrain the bulk motion of the hot gas in the flow. A crucial advantage that DN offer over other classes of systems that could harbor a hot flow is the (usually) relatively good knowledge of the mass (and therefore the size) of the accreting WD and of the inclination of the binary system. This provides an absolute scaling for the width of the emission lines, provided the line-of-sight velocities are known. Note that, except for small values of the inclination angle, broadening by the flow bulk motion is expected to dominate over thermal broadening for determining the width of the thermal X-ray lines in an ADAF.

Panels (a) and (b) of Figure 3 show the profiles of two strong lines, from iron and oxygen atoms respectively, for the inclination angle $`i=40^o`$ of SS Cyg. The predicted width of the lines is $`500`$ km s-1. Note that the oxygen line is narrower than the iron line because it is produced in cooler regions of the ADAF, where it rotates more slowly. The width of these lines is a direct probe of the sub-Keplerian rotation in the hot flow because the inclination of the system (and to some extent the region of emission) is known. In particular, if the rotation rate of the flow were Keplerian, lines would be expected to be $`5`$ times broader. Since the sub-Keplerian rotation is due to partial radial pressure support by the hot gas, and this gas is hot because of energy advection, detection of these “narrow” lines would indicate the presence of energy advection in the flow.

### 4.2 Detectability with Current X-ray Instruments

At a distance of $`166`$ pc (Harrison et al. 1999), the typical X-ray flux of SS Cyg ($`0.510`$ keV) is $`10^{11}`$ erg cm-2 s-1 or more (Mukai & Shiokawa 1993). This allows a high quality spectrum to be taken in a few tens of ks of observation with the powerful X-ray satellites Chandra and XMM-Newton. The two satellites have spectroscopic instruments which approximately cover the soft and hard X-ray energy band corresponding to the spectra shown in Fig. 2. If a radially extended ADAF is present in SS Cyg, a spectrum taken by either of them should therefore yield valuable constraints on the general structure of the hot flow (radial extent, typical densities).

An important question is whether or not the width of the predicted lines (see Fig. 3) can be resolved by current instruments. The resolving power of the HETG on Chandra at 6.9 keV is 200 (corresponding to $`1500`$ km s-1), which is of order 3 times the FWHM of the Fe XXVI line expected from the ADAF. This line cannot be resolved, but it should be possible, from a high quality spectrum, to rule out a line which is 5 times wider than predicted for the ADAF model. It should therefore be possible to test the Sub-Keplerian character of the hot flow with this instrument. Similarly, the resolving power of the HETG around 1 keV is 1000 (corresponding to $`300`$ km s-1), which is comparable to the width of the OVIII line expected from the ADAF. If such a line is detected at high enough S/N, it should be possible to tell with confidence if the hot emitting gas rotates at a substantially sub-Keplerian rate or not.

The two EPIC cameras onboard XMM-Newton have a resolving power of $`15`$ at 1 kev and $`50`$ at 6.4 keV, which is insufficient for the line diagnostics envisioned here. However, the RGS instrument (covering only the soft X-rays, from 0.35 to 2.5 keV) has a resolving power of $`300500`$ which should be sufficient to distinguish between a Keplerian and a substantially sub-Keplerian rotation rate. A possible advantage of RGS over the Chandra instruments is its larger sensitivity at low energies, allowing spectroscopy with higher time resolution.

## 5 Discussion

There is already evidence in the literature for X-ray line emission from DN in general and SS Cyg in particular. Mukai & Shiokawa (1993) find, from EXOSAT observations, that an emission line near 6.7 keV originates from many DN, including SS Cyg. Yoshida et al. (1992) find with the Ginga satellite that the X-ray emission of SS Cyg in quiescence shows an emission line around 6.68 keV (see also Done & Osborne 1997). These observations are therefore encouraging in suggesting the presence of a hot flow around the WD. Their quality, however, is insufficient to carry out the detailed spectroscopic diagnostics described above.

These diagnostics will become possible with the new powerful X-ray satellites Chandra & XMM-Newton. It should be possible, from high quality X-ray spectra, to discriminate between the case of accretion via an ADAF and the presence of a disk boundary layer (BL) in a quiescent DN. To our knowledge, there is no detailed spectral prediction for the emission originating from a compact BL (at low accretion rates), but one expects emission from a much denser gas than in the radially extended ADAF, possibly with temperatures ranging from the keV range all the way down to the WD atmospheric temperature at $`15,000`$ K or so. In first approximation, the dissipation of the rotational kinetic energy of the gas in the boundary layer leads to temperatures $`T`$ such that $`3k_BT1/2m_p(V_K^2V_{}^2)`$, where $`V_K`$ is the Keplerian speed in the disk and $`V_{}`$ is the speed at the stellar surface. If $`T`$ drops toward the WD surface, lines from cooler regions should have smaller rotational widths. For example, the model of Mahasena & Osaki (1999) does not treat the azimuthal velocity explicitly, but it predicts that most of the energy liberated in the boundary layer is transported to a dense region close to the WD surface and radiated away at intermediate temperatures. Because the radiating gas is the part of the BL closest to the stellar surface, its angular speed should be close to that of WD, and the line widths should be even smaller than in the ADAF case. However, broadening by turbulent motions and the fact that detailed and more realistic calculations show that energy is radiated radially outward and advected radially inward in BLs (e.g. Narayan & Popham 1993; Popham & Narayan 1995; Godon, Regev & Shaviv 1995; Regev & Bertout 1995) could modify this simple picture quite substantially.

A possibly important difference between an ADAF model and simple BL models (in which all the kinetic energy of the accreting gas is radiated as one-temperature Bremsstrahlung emission and the role of advection is neglected; e.g. Patterson & Raymond 1985a) is in the ionization state of the gas, as manifested in the emission line equivalent widths. The fundamental difference is that only a small fraction of the accretion luminosity in the ADAF model emerges as fairly energetic X-rays, while most of the luminosity of such a simple boundary layer does. Thus photoionization may be significant and the ionization parameter higher in the simple BL model. Elements such as oxygen will be more highly ionized, with O VII and O VIII replaced by O IX. Thus the O VII $`\lambda `$22 and O VIII $`\lambda `$19 lines should be weaker in the simple BL model than in an ADAF model. In that respect, the models presented in this study (specific to SS Cyg) have general properties which are representative of what is expected from ADAFs possibly present in other quiescent DN. More specific BL model predictions for emission line fluxes and profiles are badly needed.

The spectral models presented in this study are for specific values of the ADAF model parameters (such as $`\alpha _{\mathrm{ADAF}}=0.2`$ and $`\dot{M}=7\times 10^3\dot{M}_{\mathrm{Edd}}`$ $`1.5\times 10^{10}M_{}`$ yr-1). We find, however, that the spectral predictions would be only marginally affected, or affected in trivial ways, if these parameters were changed, so that we felt it was unnecessary to explore here the parameter space of the ADAF models. For example, line profiles are sensitive to the value of $`\alpha _{\mathrm{ADAF}}`$ only for small values of the inclination333It is only the radial component of the flow velocity that depends on $`\alpha `$, but not the rotational one. At small inclinations, turbulent motions become the main source of line broadening, together with thermal motions., while line ratios are completely insensitive to changes in $`\alpha _{\mathrm{ADAF}}`$444The run of density is independent of $`\alpha _{\mathrm{ADAF}}`$ (only the overall scaling depends on it), and the run of temperature does not directly depend on the value of $`\alpha _{\mathrm{ADAF}}`$ (it depends only on the value of $`\beta `$; see, e.g., Menou 2000).. Changes in the values of $`\dot{M}`$ do affect the overall spectrum, by rescaling the overall continuum and line emission, but without relative changes between the lines. The rate of rotation, on the other hand, depends rather sensitively on the adiabatic index $`\gamma `$ of the fluid. Note, however, that even if extreme values are assumed for this variable in an ADAF, the flow still remains sub-Keplerian ($`\mathrm{\Omega }0.5\mathrm{\Omega }_\mathrm{K}`$ when $`\gamma 4/3`$; Narayan & Yi 1995) and our conclusions are mainly unchanged.

A possible ingredient of the X-ray spectrum of a quiescent DN that has been omitted from our calculations is the iron $`K_\alpha `$ fluorescence line at $`6.4`$ keV from the X-ray irradiated disk. This line can appear on top of the Bremsstrahlung continuum emission, depending on the radial extent of the truncated disk. The width of the line, related to the Keplerian rotation rate in the inner regions of the disk and the system inclination is $``$ a few $`10`$ eV at most (Menou 2000). This is still narrow enough that the line should be easily separable, in a good quality spectrum, from the other iron lines expected at $`6.7`$ and $`6.9`$ keV (from the hot gas).

Finally, it is important to notice that there are possible complications in interpreting narrow X-ray emission lines in terms of energy advection in a hot flow. One possible source of confusion comes from lines emitted in the close vicinity of the WD (for instance in a disk BL scenario), which is also likely to be rotating at a sub-Keplerian rate (see Sion 1999 for existing constraints on the rotation rate of WDs in DN). The sub-Keplerian rotation of the emitting gas can therefore be associated to the property of energy advection of the hot flow only if it can be shown that it is not directly related to the slow rotation of the WD. This may require clear evidence that the hot flow is indeed radially extended. We also note that Medvedev & Narayan (2000) recently presented solutions for accretion onto a compact object in which the flow can be radially extended and have a substantially sub-Keplerian rotation rate, without advecting any energy.555The rotation rate is directly related to that of the central object in these solutions. If such a hot flow were present around WDs in quiescent DN, it could also be at the origin of narrow X-ray emission lines. The interpretation of narrow lines in terms of energy advection will therefore require some knowledge of other properties of the flow, such as its density profile, which differs in the ADAF solutions and the solutions of Medvedev & Narayan (2000). The knowledge of the WD rotation rate in some DN and the possibility of actually resolving some of the low energy X-ray lines predicted for ADAFs (§4.2) should provide critical observational tests of these two accretion scenarios.

## 6 Conclusion

We have constructed detailed spectral models of an ADAF around a WD, with model parameters appropriate for the Dwarf Nova SS Cygni in quiescence. The ADAF is a source of X-ray Bremsstrahlung continuum emission, which is characteristic of its density, temperature and radial extent.

Strong X-ray lines are also emitted by highly ionized atoms of Fe, O, Si and other heavy elements in the flow. The pattern of emission lines is a robust probe of the density, temperature and elemental abundances in the hot flow, which can therefore be constrained from high quality X-ray spectra.

In addition, the width of the X-ray emission lines can be directly related to the bulk motion of the gas in the ADAF by using the a priori knowledge of the mass of the WD and the inclination to the line of sight of the system. Current X-ray instruments have the capabilities to distinguish between lines emitted from a flow rotating at or substantially below the Keplerian rate. The property of energy advection of the flow, which is at the origin of the sub-Keplerian rotation in an ADAF, can therefore in principle be tested accurately by the observations.

## Acknowledgments

We are grateful to Jonathan McDowell for many enlightening discussions on the observational capabilities of Chandra. Support for this work was provided by NASA through Chandra Postdoctoral Fellowship grant number PF9-10006 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-39073. Partial support was provided by NASA Grant NAG-528 to the Smithsonian Astrophysical Observatory.


format: text

text:
# 3+1 dimensional Yang-Mills theory as a local theory of evolution of metrics on 3 manifolds.

## A Appendix

To evaluate $`\frac{\delta \theta ^b[e]}{\delta e_i^a}`$, note that

$$e_i^a+\delta e_i^a=(e_{ij}+\delta e_{ij})exp(i(\theta [e]+\delta \theta [e])^b\frac{T^b}{2}).$$
(34)
This gives

$$\delta e_i^aO^{ka}=\delta e_{ik}+ϵ^{dae}e_{ij}O^{jd}O^{ka}\overline{\delta \theta }^e$$
(35)
where $`\overline{\delta \theta }^e=M_e^f[e]\delta \theta ^e`$ and the matrix $`M_e^f`$ is given by

$$O^T(\theta )O(\theta +\delta \theta )1+T^aM_b^a(\theta )\delta \theta ^b.$$
(36)
Taking the antisymmetric part of (35) we can solve for $`\frac{\delta \theta ^b[e]}{\delta e_i^a}`$. We get

$$\frac{\delta \theta ^b[e]}{\delta e_i^a}=ϵ_{ijk}M^1[\theta ]_c^b((e\mathrm{𝟏}tre)^1)_{kl}O^{ja}O^{lc}.$$
(37)
When two or all eigenvalues of the symmetric matrix $`e_{ij}`$ are degenerate, so are those of $`(e\mathrm{𝟏}tre)`$. In this case, the variables $`\theta [e]`$ equation (34) are ill defined and more care is required to define the new variables.


text:
# Minimum Entropy Approach to Word Segmentation Problems

## I Introduction.

The problem addressed in this paper is rather elementary in statistics. It is best described as the following: suppose one who knows nothing about English language was given a sequence of English letters, which was actually obtained by taking off all the interwords delimiters among a sample of English text, how could he recover the words of the text by choosing to insert spaces between adjacent letters? Note that the only thing he can consult is the statistical properties of the sequence?

Any two adjacent letters can be chosen to belong to the same word (keep adjacent) as well as belong to separate words (be separated by space). Suppose the sequence length is $`N`$. Any choice on the connectivity between $`N1`$ pairs of adjacent letters is called a segmentation. There are a total of $`2^{N1}`$ possible segmentations. The word segmentation problem is to find ways to distinguish the correct segmentation – in the sense that adjacent letters in the original text keep adjacent while letters separated by spaces and/or punctuation marks in the original text are separated by spaces in the segmentation – from others.

Although the problem seems toy-like, its fundamental importance for statistical linguistics is evident. We study on it, however, also for practical purposes. Noncoding sequences in the genomes of species play essential rule on the regulation of gene expression and function . However the development of computational methods for extracting regulatory elements is far behand DNA sequencing and gene finding . One reason is the lack of efficient way to discriminate large amount of sequence signals in noncoding DNA sequences. Through linguistic study it has been shown that noncoding sequences in eukaryotic genomes are structurally much similar to natural and artificial language . Thus many may expect to “read” the noncoding sequences as a “text”. Actually, efforts have been given to build a dictionary for genomes . Li et al.  showed the connection between regulatory elements prediction and word segmentation in noncoding DNA sequences of eukaryote genomes. We expect that progress on word segmentation problem may help to deepen our knowledge on noncoding regions of eukaryote genomes. Besides, word segmentation is an important issue for Asian languages (e.g., Chinese and Japanese) processing , because they lack interword delimiters.

## II Segmentation entropy and its connection to word segmentation problem.

To tackle word segmentation problem, we first consider a problem under constraints, so that one important concept – segmentation entropy – can be introduced. The constraints will be released at the end of this paper. Suppose we have known that there are $`n_l`$ words of length $`l`$ $`(l=1,2,\mathrm{})`$ in the original text. Obviously,

$$\underset{l}{}n_ll=N.$$
(1)
Under these constraints – Words Length Constraints WLC – there are totally

$$\frac{(_ln_l)!}{_ln_l!}$$
(2)
segmentations. For example, for the following story, there are totally $`3.12e144`$ segmentations, while the number under WLC is about $`1.33e97`$.

> The Fox and the Grapes
>
> Once upon a time there was a fox strolling through the woods. He came upon a grape orchard. There he found a bunch of beautiful grapes hanging from a high branch.
>
> “Boy those sure would be tasty,” he thought to himself. He backed up and took a running start, and jumped. He did not get high enough.
>
> He went back to his starting spot and tried again. He almost got high enough this time, but not quite.
>
> He tried and tried, again and again, but just couldn’t get high enough to grab the grapes.
>
> Finally, he gave up.
>
> As he walked away, he put his nose in the air and said: “I am sure those grapes are sour.”

Following least effort principle , it is appreciable in natural languages to combine existing words to express different meaning. Shannon  pointed out the importance of redundancy in natural languages long ago: generally speaking, nearly half of the letters in a sample of English text can be deleted while someone else can still restore them. These properties of natural language ensure the sequence obtained by taking off interword delimiters from a certain text being highly nonrandom and showing determinant and regular characteristics. It is expected that the correct segmentation reveals these characteristics as strongly as possible. From information point of view, this means that, if a form of information entropy can be properly defined on each segmentation, the entropy of the correct segmentation will be the smallest.

Interestingly, a maximum likelihood approach leads to the same proposal and automatically gives the definition of the entropy. Given one sequence of length $`N`$, we expect to find a likelihood function which reaches its maximum on the correct segmentation. For a concrete segmentation, we assign a probability to each word in it

$$w_ip_i,i=1\mathrm{}M$$
(3)
with

$$\underset{i=1}{\overset{M}{}}p_i=1.$$
(4)
The likelihood function is written as

$$Z_s=\underset{i=1}{\overset{M}{}}p_{i}^{}{}_{}{}^{m_il_i}$$
(5)
where $`m_i`$ is the number of word $`w_i`$ in the segmentation, and $`l_i`$ is the length of the word.

By maximizing the likelihood function subjected to eq.(4) we obtain

$$p_i=\frac{m_il_i}{N}.$$
(6)
Thus the maximum likelihood for the segmentation is

$$Z_s=\underset{i=1}{\overset{M}{}}(\frac{m_il_i}{N})^{m_il_i}.$$
(7)
The segmentation with maximum likelihood is just the one minimizing

$$S=\frac{lnZ_s}{N}=\underset{i=1}{\overset{M}{}}\frac{m_il_i}{N}ln(\frac{m_il_i}{N}).$$
(8)
This function has the form of entropy  and will be called Segmentation Entropy (SE).

Starting from a maximum likelihood approach, we now come to the suggestion to minimize the segmentation entropy. This is in contrast to commonplace. Maximizing likelihood leads to maximizing certain entropy in some cases . As a general principle for investigating statistical problems, maximum entropy method has been successfully applied in a variety of fields . We propose that, instead of applying maximum entropy principle, one may choose to minimize certain entropy (minimum entropy principle) in some problems. This seems attractive especially in the era of bioinformatics when most of the problems are to reveal regularity in large amount of seemingly random sequences.

Because the present is a statistical method, the text under study needs to be not too short. For example, when we tried to find the segmentation with the smallest segmentation entropy for the saying

> God is nowhere as much as he is in the soul… and the soul means the world

(By Meister Eckhart, 14-century Dominican priest, Preacher, and Theologian), it was found that, among a total of $`343062720`$ segmentations under WLC, there are 15 segmentations whose SE is $`2.3684`$, smaller than 2.3802 of the correct one. One example is

> god isnow he rea smuchas he is int he soul andt he soul meanst he world,

in which the five “he” and two “soul” are revealed.

Unfortunately, present computational power does not permit to exhaustively study even a text as short as “the Fox and the Grapes”, the number of permitted segmentations for which is $`1.33E+97`$ under WLC. We choose to see the relevance of the concept of segmentation entropy in some special ways. The study focuses on “The Fox and the Grapes”.

To change a segmentation slightly, one way is to choose two adjacent words along the sequences randomly and then exchange their length. This way the original two words may change to different words. This procedure can be repeated on the resulting segmentations. The change does not violate the WLC. Because of the large number of possible choices in each step, the segmentation is expected to become increasingly dissimilar to the original one. Starting from the correct segmentation of “The Fox and the Grapes”, we expect to see the evolution of SE by changing the segmentation this way. Figure 1 shows that SE increase drastically in the first 500 steps, and then reaches and fluctuates around certain equilibrium value. Compared with the gap between the equilibrium value and the original SE, the fluctuation is minor. This shows that, at least locally, the correct segmentation is at the minimum of SE. Actually, we have traced a trajectory of evolution up to $`10^{10}`$ steps. No segmentation with SE smaller than the correct one was observed. This implies that SE of the correct segmentation is also globally minimal.

The distribution of segmentation entropy may give further insight to the atypicality of the correct SE. We randomly sampled $`10^{10}`$ segmentations in the following way: while keeping the WLC, the length of each words in the segmentation is assigned randomly. The distribution of SE is shown in Fig. 2. The minimal SE we sampled is 4.5298, still much higher than 4.097 of the correct segmentation (see Fig. 1). It is interesting to observe that the distribution shows fractal characteristics. The fractal-like distribution presents also for other text, even for random sequence (Fig. 3). The fractal-like feature is determined by the WLC and the statistical structure of the sequence under study. In Fig. 3 we compared the distribution of SE of two sequences (under the same WLC), the original sequence of “The Fox and the Grapes” and a random sequence obtained by randomizing the order of letters in the text. The result is in accordance with the fact that the original sequence is in a much more ordered state, manifesting that segmentation entropy captures the statistical structure of the sequences successfully.

There is one way to estimate the number of segmentations the SE of which is 4.097, the value for the correct segmentation. See Fig. 4 in which the distribution of SE in Fig. 2 are shown in logrithmic scale here. The left edge of the distribution fall on a line. The edge can be fitted by $`e^{(165x750.42)}.`$ The number of segmentations with SE x among the totally $`1.33e97`$ possible segmentations under WLC is:

$$c(x)=\frac{1.33e^{97}}{9\times 10^9}e^{(165x750.42)}.$$
(9)
We obtained $`c(4.097)=0.96`$. From the distribution of SE shown in Fig. 3(a) we obtained the same value of $`c(4.097)`$. The estimation support the idea that segmentation entropy of correct segmentation is unique.

We now consider how to release the WLC. Unfortunately, searching the segmentation with the smallest SE among all the possible is sure to fail to find the correct one. For example, SE of the segmentation in which the whole sequence is considered as one word (single-word segmentation) is 0, the smallest possible SE. Also, the segmentation in which each letter is viewed as a separate word ($`N`$-word segmentation) has a considerably small SE (2.8655 for “The fox and the grapes”). These are called side attraction effects. These examples show that smaller SE does not necessarily means better segmentation when we compare the SEs of segmentations under different WLC (here WLC refers to any partition of numbers of words of various length satisfying eq.(1), not necessarily the same as the original text.) The bias induced by different WLC must be taken off. In order to do so, we suggest to use

$$R_S=\frac{S}{S_0}$$
(10)
instead of $`S`$. Here $`S_0`$ is the average SE under the same WLC of a sequence obtained by randomizing the order of letters in the original text. $`S_0`$ plays the role of chemical potential for a thermodynamic system . $`R_S`$ for the single word and $`N`$-word segmentations are 1, the largest possible value. By searching segmentation with the smallest $`R_S`$, it is expected to find meaningful segmentation. For examples, for the segmentation

> god isnow he rea smuchas he is int he soul andt he soul meanst he world,

which has already been shown above, $`R_S`$ is 0.8601; while

> god is now he re as much as he is int he soul an dt he soul me an st he world

is a better – actually one of the best – segmentation according to $`R_S`$ ($`R_S=0.8259`$). Intuitively this is reasonable, because in this second segmentation, more repeated “words” – two copies of “is”, “as” and “an” – are revealed. Another segmentation

> god is now he re as much as he is in thesoul an d thesoul me an st he world,

which differs from the second segmentation by revealing the two “thesoul”, has a moderately small $`R_S`$: 0.8481. Comparison shows that the five repeats of “he” is the most preferred part in good segmentations.

## III Concluding remarks.

In statistical linguistics many efforts are given on signal extracting and statistical inference. Our method, however, is new on at least two points. First, there is neither assumption on distribution  nor demand for training sets, lexical or grammatical knowledge . This feature is important for studying biological sequences, because present knowledge on the “language” (DNA) of life is still lack. Second, instead of extracting a limit number of signals, we try to “read” the sequence exactly as a “text”. A text includes more than words: it also includes the organization of words. The results of segmentation form a basis for many further elaborations.

Principally, the concept of segmentation entropy can be applied to study the noncoding DNA sequences of eukaryote genomes. It is expected that the study may gives more than some meaningful “words” or regulatory elements. Possible applications are not confined to studying noncoding DNA sequences of course. Segmentation entropy can be used to find patterns in any symbolic sequences. However, the application of segmentation entropy is restricted by the difficulty to find the segmentation with the smallest $`R_s`$ from the vast amount possible ones. We are now developing algorithm that can be used for regulatory binding sites prediction. in the algorithm the principle of minimun entropy will be incorporated in.

## ACKNOWLEDGMENTS

I thanks Professor Bai-lin Hao who helps to make the computing possible. I also thanks Professor Wei-mou Zheng and Professor Bai-lin Hao for stimulating discussions. Mr. Xiong Zhang carefully read the manuscript. The work was supported partly by National Science Fundation.


text:
# Tidally–induced angular momentum transport in disks

## 1. Angular momentum exchange between the disk rotation and the orbital motion

In a binary system where at least one of the stars is surrounded by a circumstellar disk, tidal waves excited by the companion propagate into the disk. If the disk and the orbital plane are coplanar, these waves are called density waves. In a noncoplanar system, both density and bending waves are excited. They are respectively of even and odd symmetry with respect to reflection in the disk midplane. The pattern speed $`\mathrm{\Omega }_P`$ with which the tidally excited pattern rotates is $`\omega `$ and $`2\omega `$ for density and bending waves, respectively, where $`\omega `$ if the binary angular speed. The disk is truncated by tidal effects in such a way that its radius is not greater than about one–third of the separation of the system (Papaloizou & Pringle 1977, Paczyński 1977, Larwood et al. 1996 for the non coplanar case). Therefore, $`\mathrm{\Omega }_P`$ is smaller than the angular speed of the gas in the disk, and the tidal pattern carries negative angular momentum (or, in other words, the torque exerted on the disk by the companion is negative). Through dissipation of the waves, the disk then loses angular momentum and disk material flows inwards, whereas the companion star, which excites the waves, gains this angular momentum.

Since there is no corotation resonance in the disk, secular exchange of angular momentum between the disk rotation and the orbital motion occurs only if the waves dissipate, either in the bulk of the disk or at its boundaries (Goldreich & Nicholson 1989). In a laminar disk, dissipation of tidal waves may arise through shocks. A shock front forms when the group velocity of the wave relative to that of the fluid (or, equivalently, the perturbed velocity) is supersonic. Shocks are very dissipative so the wave amplitude cannot rise much above the level where the front first forms. In other words, the wave is decelerated at the shock front in such a way as to restore marginally sonic wave motion. Conservation of wave action may tend to cause the amplitude of the (shock–)wave to increase again as it propagates further in, but this effect is balanced by shock dissipation maintaining the amplitude at the marginal level.

Tidally–induced angular momentum transport has been considered as an alternative to turbulent transport (Shu 1976, Sawada et al. 1986, Spruit et al. 1987, Larson 1989, see also Larson in this volume). It is viewed as particularly attractive in disks where the ionization level is too low for the magnetorotational instability to operate (Balbus & Hawley 1998). However, while the presence of spiral waves have been inferred in the accretion disk of the dwarf nova IP Peg (Steeghs et al. 1997), there are no indications that they are associated with significant angular momentum transport. Below we consider successively the case of density and bending waves.

## 2. Angular momentum transport by density waves

There are two main obstacles to the propagation of density waves down to small radii. First, if the disk is vertically stratified, the wave front tends to be tilted upwards so that the wave action migrates towards the surface of the disk and into any atmosphere it possesses where it can take on high amplitude and be dissipated (Lin et al. 1990a, 1990b). Only under the artificial assumptions of a strict polytropic edge and no dissipation can it be channeled into and remain in a very narrow surface waveguide (Ogilvie & Lubow 1999). Migration of wave action towards the surface (or wave refraction) is more effective for high azimuthal mode number $`m`$, but is still efficient for the two–armed spirals which are predominantly excited in binaries. This is because even though these waves have a larger wavelength in the linear regime, if they become nonlinear their profile necessarily distort and develop short wavelength components for which refraction might be important.

The second obstacle to long range wave propagation is a low temperature or, equivalently, a high Mach number (e.g., Spruit 1987, Savonije et al. 1994, Godon et al. 1998, Armitage & Murray 1998, Blondin 2000). This is because the characteristic wavelength of the excited waves decreases with increasing Mach number, so that it becomes very small compared with the scale associated with the forcing potential (Lin & Papaloizou 1993). The torque, which is obtained by integrating over the volume of the disk the perturbed mass density times the tidal force, is then very small.

There have been a number of 3D numerical calculations of tidal shock waves (see Yukawa et al. 1997 and references therein, Haraguchi et al. 1999) but the loss of disk angular momentum has not been computed in these calculations. Spiral shocks were seen in some of these simulations, but they were much less distinct than in 2D. Even in 2D, where refraction is absent and therefore tidal effects are overestimated, the pattern observed in IP Peg can be reproduced only in the outer disk during outburst, when the enhanced viscosity pushes the disk edge into a region of strong gravitational perturbations from the secondary (Armitage & Murray 1998), or for unrealistically hot disks (Godon et al. 1998). Note that observations themselves only show spirality in the outer disk in IP Peg. In 2D, calculations by Savonije et al. (1994) and Blondin (2000) suggest that wave–driven accretion onto the central star occurs in disks where the Mach number is smaller than about 10, whereas it is inefficient when the Mach number is larger than about 20. Since the Mach number in protostellar disks is thought to be larger than 10 (the disk aspect ratio is around 0.05–0.1), tidally–induced transport is probably not significant in the disk inner parts.

## 3. Angular momentum transport by bending waves

Bending waves are more efficient at transporting angular momentum in a disk than density waves, because they have a longer wavelength (Papaloizou & Lin 1995). For the same reason they can also propagate down to smaller radii.

Papaloizou & Terquem (1995) calculated the $`m=1`$ bending wave response of an inviscid disk with the rotation axis misaligned with a binary companion’s orbital rotation axis. They assumed that the waves were dissipated by nonlinear interaction with the background flow before reaching the disk inner edge, so that all their angular momentum was deposited into the disk. They found that $`m=1`$ bending waves can lead to the accretion of the disk on a timescale in excess of a few times $`10^7`$ years.

Terquem (1998) calculated the tidal torque in a viscous disk (where the waves are viscously damped) and found that it can be comparable to the torque communicated internally by horizontal viscous stress acting on the background flow when the perturbed velocities in the disk are on the order of the sound speed. The tidal torque can exceed the horizontal viscous torque only if the viscous stress tensor is anisotropic with the parameter $`\alpha `$ which couples to the vertical shear being larger than that coupled to the horizontal shear. When the perturbed velocities become supersonic, shocks reduce the amplitude of the perturbation such that the disk moves back to a state where these velocities are marginally sonic (Nelson & Papaloizou 1998). When shocks occur, the tidal torque exerted on the disk may become larger than the horizontal viscous torque. Terquem (1998) also found that if the waves are reflected at the center, resonances occur when the frequency of the tidal waves is equal to that of some free normal global bending mode of the disk. If such resonances exist, tidal interactions may then be important even when the binary separation is large. However, it is unlikely that in a realistic accretion disk waves can be reflected at the disk inner edge. Therefore, in a viscous disk, it is unlikely that transport of angular momentum is increased by more than a factor two or so by tidal effects. It was also found that if $`\alpha `$ is larger than about $`10^2`$, bending waves are damped before they can reach the disk inner parts.

## 4. Conclusion

Tidal effects in pre–main sequence binary systems are important for truncating protostellar disks and for determining their size. However, once the disk is truncated, the calculations reviewed above suggest that tidally–induced angular momentum transport does probably not play a dominant role. Where the disk is ionized enough so that the magnetorotational instability can develop, transport of angular momentum is most probably dominated by magnetic turbulence. If the disk is laminar, the amount of transport provided by tidal waves is unlikely to be large enough to account for the dissipation of the disk on a timescale on the order of a few million years. In addition, tidal effects tend to be exerted mainly at the disk outer edge, where the perturbation is the strongest. Strong tidal effects at the outer edge of the disk would allow mass from the outer region to retract inwards, subsequently weakening the tides at the edge. Some pile up of mass may result at smaller radii, but whether that would affect significantly the whole disk would depend on the disk mass.

### Acknowledgments.

I thank S. Balbus and J. Papaloizou for useful comments on an early draft of this paper.

## References

Armitage, P. J., & Murray, J. R. 1998, MNRAS, 297, L81

Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1

Blondin, J. M. 1999, preprint (astro–ph/9909181)

Godon, P., Livio, M., & Lubow, S. 1998, MNRAS, 295, L11

Goldreich, P., & Nicholson P. D. 1989, ApJ, 342, 1075

Haraguchi, K., Boffin, H. M. J., & Matsuda, T. 1999, in Star Formation 1999, ed. T. Nakamoto (Nobeyama Radio Observatory), p. 241

Larson, R. B. 1989, in The Formation and Evolution of Planetary Systems, eds. H. A. Weaver & L. Danly (Cambridge: CUP), p. 31

Larwood, J. D., Nelson, R. P., Papaloizou, J. C. B., & Terquem, C. 1996, MNRAS, 282, 597

Lin, D. N. C., & Papaloizou, J. C. B. 1993, in Protostars and Planets III, ed. E. H. Levy & J. Lunine (Tucson: Univ. Arizona Press), p. 749

Lin, D. N. C., Papaloizou, J. C. B., & Savonije, G. J. 1990a, ApJ, 364, 326

Lin, D. N. C., Papaloizou, J. C. B., & Savonije, G. J. 1990b, ApJ, 365, 748

Nelson, R. P., Papaloizou, J. C. B., 1999, MNRAS, 309, 929

Ogilvie, G. I., & Lubow, S. H. 1999, ApJ, 515, 767

Paczyński, B. 1977, ApJ, 216, 822

Papaloizou, J. C. B., & Lin, D. N. C. 1995, ApJ, 438, 841

Papaloizou, J. C. B., & Pringle, J. E. 1977, MNRAS, 181, 441

Papaloizou, J. C. B., & Terquem, C. 1995, MNRAS, 274, 987

Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13

Sawada, K., Matsuda, T., & Hachisu, I. 1986, MNRAS, 219, 75

Shu, F. H. 1976, in Structure and Evolution of Close Binary Systems, IAU Symp. 73, eds. P. Eggleton, S. Mitton, & J. Whelan (Reidel: Dordrecht), p. 253

Spruit, H. C. 1987, A&A, 184, 173

Spruit, H. C., Matsuda, T., Inoue, M., & Sawada, K. 1987, MNRAS, 229, 517

Steeghs, D., Harlaftis, E. T., & Horne, K. 1997, MNRAS, 290, L28

Terquem, C. E. J. M. L. J. 1998, ApJ, 509, 819

Yukawa, H., Boffin, H. M. J., & Matsuda, T. 1997, MNRAS, 292, 321


text:
# Diffraction-limited near-IR imaging at Keck reveals asymmetric, time-variable nebula around carbon star CIT 6

## 1 Introduction

As evidenced by its primary designation, CIT 6 (= RW LMi = GL 1403 = IRC +30219 = IRAS 10131+3049) was discovered as a bright infrared source without a known visible counterpart at the time of the Caltech infrared survey (Ulrich et al. 1966). The identification of CIT 6 as a carbon star is based on its near-infrared CN bands (Wisniewski et al. 1967; Gaustad et al. 1969) and its thick carbon-rich dust shell, responsible for the heavy visible obscuration (e.g., Strecker & Ney 1974). Monitoring of the infrared variability, first recognized by Ulrich et al. (1966), has led to an established period of $``$628 days for this long-period variable (Alksnis & Khozov 1975; Alksnis 1995; Taranova & Shenavrin 1999). Based on the pulsational period, the distance to CIT 6 has been estimated by Cohen & Hitchon (1996) as 400$`\pm `$50 pc from application of the period-luminosity relations established for O-rich Miras (LMC and galactic) and C-rich Miras (LMC only).

Both the visible and infrared flux from CIT 6 is highly polarized, indicating an asymmetric distribution of circumstellar material (Kruszewski 1968; Dyck et al. 1971). In addition, the fractional polarization is variable (Kruszewski 1971; Dyck et al. 1971; Kruszewski 1973; Alksnis & Khozov 1975) and the polarization angle is currently strongly wavelength-dependent, changing $``$90° between the blue and the red (Kruszewski 1973; Cohen & Schmidt 1982; Trammell et al. 1994). This distinctive behavior is likely linked to the excess UV/blue flux and unusual emission lines which have been observed, perhaps due to the presence of an unseen (hot) companion (Alksnis et al. 1988) or circumstellar shocks (Cohen 1980; Cohen & Schmidt 1982).

Although no previous high resolution imaging of this source has been reported, the characteristic sizes of the near- and mid-infrared emission regions have been measured. Dyck et al. (1984) used speckle interferometry to estimate the 2.2µm, 3.5µm, and 4.8µm sizes to be $``$70 mas, 100 mas, and 144 mas respectively (Gaussian FWHM). In addition, Lipman (1998) reported a 11µm FWHM of 150 mas using the Infrared Spatial Interferometer.

Recently developed observing techniques, from both the ground and space, allow the asymmetries of the dust envelope to be observed directly. In this paper we report multi-epoch, diffraction-limited observations of the circumstellar dust shell at 1.6µm, 2.2µm, and 3.1µm from the Keck-I telescope. Using aperture masking interferometry, we attain similar spatial resolution in the infrared to that of the Hubble Space Telescope in the visible, and the brightness distributions at these different wavelengths have been registered and compared in order to probe the circumstellar environment of this enigmatic system.

## 2 Observations

### 2.1 Infrared Interferometry with Keck-I

Aperture masking interferometry was performed by placing aluminum masks in front of the Keck-I infrared secondary mirror. This technique converts the 10-m primary mirror into an interferometric array, allowing the Fourier amplitudes and closure phases for a range of baselines to be recovered with minimal “redundancy” noise (e.g., Baldwin et al. 1986). In this experiment, both non-redundant and doubly-redundant pupil configurations were employed. The Maximum Entropy Method (MEM) (Gull & Skilling 1983; Sivia 1987) was used to reconstruct diffraction-limited images from the interferometric data. MEM does produce structure in the maps beyond the formal diffraction limit, an extrapolation based on the positivity and finite-support image constraints, as well as the “maximum entropy” regularization scheme. In some cases, this structure can be spurious and result in well-known MEM artifacts, such as ringing of point sources embedded in nebulosity (Narayan & Nityananda 1986). For CIT 6, strong (non-zero) closure phase signals suppress these artifacts and result in a high fidelity map. In order to check the reliability of the reconstructions, MEM results were compared to the results of the CLEAN reconstruction algorithm (Högbom 1974). Results were always consistent, with smoothed MEM maps being comparable to the the results of CLEAN, which does not generate similar super-resolution and hence has less fine structure in its images. Further engineering and performance details may be found in Monnier (1999) and Tuthill et al. (2000b).

CIT 6 was observed at 6 different epochs between 1997 January and 2000 January, at the Keck-I telescope using the Near-Infrared Camera (Matthews & Soifer 1994; Matthews et al. 1996) and with a variety of aperture masks. Important observing information can be found in Table 1, including spectral filter characteristics, the number of speckle frames, and the unresolved stars used for calibrating the atmosphere plus telescope transfer function. In all cases, the integration time of each frame was 0.137 s.

Figure 1 shows MEM reconstructed images of CIT 6 at 2.2µm and 3.1µm from all epochs. In cases where multiple observations exist at a given epoch, images were reconstructed from the separate datasets and then averaged together, weighted by the dynamic range of each image. In addition, the last known direction of the polarized component of the electric field in the K band has been included in the first panel for easy reference (Dyck et al. 1971). At both colors the dominant component, appearing at the map center, appears as a bright feature, elongated to the West-Northwest at a position angle of about $`290^{}`$. A second component, visible but faint in the 2.2 $`\mu `$m images, and more prominent at 3.1 $`\mu `$m, appears separated by $``$60 milliarcseconds (mas) at a position angle (PA) of $``$215°from the map center. These components will be referred to as the northern and southern components respectively. These data will be discussed further in the next section.

While CIT 6 was observed near 1.6µm, low signal-to-noise ratio prevented accurate image reconstructions. However, from the visibility data, it is clear that the dominant component is compact and elongated along PA $``$286°, similar to that seen at 2.2 µm and 3.1µm. Fitting a two-dimensional Gaussian to the visibility data yielded a major axis of $``$24 mas and minor axis of $``$16 mas. This simple model does not adequately fit the data ($`\chi ^2`$ per degree of freedom $``$5), and can only be viewed as a crude description of the actual brightness distribution.

### 2.2 Hubble Space Telescope images

CIT 6 has also recently been imaged with the Hubble Space Telescope (HST) using both WFPC2 and NICMOS instruments. The WFPC2 observations were carried out on 1996 April 2, and the reduced images discussed here were provided by S. Trammell in advance of publication (Trammell et al. 2000). The plate scale of the original images was 46 mas per pixel. Figure 2 shows the logarithmic intensity of CIT 6 using filters F439W ($`\lambda `$=0.429µm, $`\mathrm{\Delta }\lambda `$=0.046µm), F555W ($`\lambda `$=0.525µm, $`\mathrm{\Delta }\lambda `$=0.122µm), and F675W ($`\lambda `$=0.674µm, $`\mathrm{\Delta }\lambda `$=0.089µm). See the WFPC2 Instrument Handbook for more information on the filter bandpasses ($`\mathrm{\Delta }\lambda `$ is roughly the FWHM). The polarization of the electric field in each of these bandpasses has been included in the figure, both Cohen & Schmidt (1982) and Trammell et al. (1994) results. Because of strong wavelength-dependence of the polarization direction within the large bandpass of the F555W filter, two polarization vectors for each epoch have been plotted to represent the range of observed angles.

The reddest image (F675W) consists of two components, or lobes, of similar brightness, separated by $``$190 mas along PA 190°. Moving towards shorter wavelengths, the F555W image still shows evidence for two lobes, but the northern component is considerably fainter with respect to the southern lobe. Lastly, the F439 image shows only the southern component, exhibiting a peculiar comet-like asymmetry at a position angle roughly orthogonal to the two lobes seen in red. Rayleigh scattering along this PA would explain the different angle of polarization for the blue light, compared to the red emission (e.g., Cohen & Schmidt 1982). Detailed analyses of these images can be found in Trammell et al. (2000).

CIT 6 was observed multiple times by NICMOS before the cryogen expired, and the data from 1998 February 23 and 1998 May 5 were used here. The calibrated data were retrieved from the HST Data Archive, including the observation log files. The images themselves will not be discussed in this paper, but rather the data will be used for extracting relative astrometry. This has allowed the Keck infrared images to be registered with respect to the WFPC2 images as is discussed fully in the next section.

### 2.3 Image Registration

Even though the Keck imagery contains no positional information and HST absolute astrometry is only good to $``$1″, it has been possible to link the infrared and visible reference frames to about 35 milliarcseconds. This is because CIT 6 has a companion star located only $``$10.4″ away (Bastian et al. 1988, hereafter referred to as the reference star), allowing precise relative astrometry. In this section we describe the methodology for extracting precise relative astrometry from both WFPC2 and NICMOS images.

#### 2.3.1 WFPC2

Deriving the position of CIT 6 relative to the reference star is straightforward because both appear on the PC chip in all of the WFPC2 images, which has a nominal plate scale of 45.5 mas/pixel. However, for the images taken with filters F439W and F555W, the reference star is slightly saturated, reducing the precision of the centroiding process. When not saturated, the centroid was defined by fitting a pixelized, two-dimension Gaussian profile to a 3x3 pixel box centered on the peak pixel. When saturated, the centroid was defined by fitting a Gaussian to the non-saturated pixels contained in a 7x7 pixel box. Tests were performed on unsaturated point spread functions and differences in the estimated centroids resulting from these two methods were always less than 6 mas. The IRAF-STSDAS task metric was translated into IDL and used to correct for the geometric field distortion (Gilmozzi et al. 1995), however corrections were always $`\stackrel{<}{_{}}`$0.3 pixels ($``$13.7 mas). For F439W only the south component was visible, but both components of the CIT 6 nebula were measured relative to the reference source for F555W and F675W images. Although the HST Data Handbook claims that this process should be accurate to 5 mas for relative astrometry on the same chip, we have assigned more pessimistic error estimates because the centroid was not as well defined for the distorted isophotes of the extended emission, and because of saturation in the core of the reference star, most pronounced in F555W. Relative astrometry with estimates of uncertainty are included in Table 2 and the results are plotted in Figure 3.

#### 2.3.2 NICMOS

NICMOS observations of CIT 6 have been carried out (by other workers, P.I. Gary Schmidt), and reside on the HST Data Archive. Although the spatial resolution of HST is four times worse than the Keck-I Telescope, these data are useful in the context of this paper for linking the position of the infrared peak to the visible WFPC2 images. Geometric field distortion was corrected according to Cox et al. (1997), while the appropriate platescale was determined from the online calibration file (http://www.stsci.edu/instruments/nicmos/NICMOS\_status/platescale.dat). When the images were not saturated, the centroid was determined as above, using a pixelized Gaussian to fit to a 3x3 pixel square surrounding the peak pixel. However, two practical problems complicated the analysis.

Firstly, the field-of-view of the NIC1 chip (used for J and H band observations) is only 11″$`\times `$11″; not large enough to image both CIT 6 and the reference star without placing CIT 6 on the very edge of the chip. However, high precision relative astrometry can be maintained by executing a “small angle maneuver.” When the telescope pointing is changed using this procedure, the angular offset between images should be known to $``$10 mas, since the same two guide stars are being used (G. Schneider, private communication). By using the RA\_AVG, DEC\_AVG, and ROLL\_AVG keywords in the JIF observing log file, the relative offsets between nearby pointings can be determined. We applied this procedure to each exposure of an observation to measure the relative separation between CIT 6 and the reference star in filters F110W and F160W for the observations of 1998 May 5. This method was tested using F220M images taken before and after a small angle maneuver, both of which contain the reference source. Two epochs of such data exist (see below) and the residual astrometric error was found to be only $`\stackrel{<}{_{}}`$12 mas for these test cases.

A second complication was that the F160W and F220M images of CIT 6 were heavily saturated. In these cases, the location of the CIT 6 centroid was estimated by measuring the intersection point of the diffraction spikes. We tested the accuracy of this method by applying it to the unsaturated F110W image and found an offset of only $`\stackrel{<}{_{}}`$13 mas with respect to the centroid determined by a Gaussian fit. However, since the diffraction pattern is wavelength-dependent, this does not assure such accuracy for longer wavelengths. We have assigned a conservative positional uncertainty of 0.5 pixels (22 mas for NIC1, 35 mas for NIC2) to measurements made using this method.

The NIC2 chip has a larger field of view (19.2″$`\times `$19.2″) which allowed both reference star and CIT 6 to appear on the same chip for the F220M data. Hence, the relative astrometry in this waveband does not suffer from possible offset errors associated with the small angle maneuver, but only those due to saturation of CIT 6 and uncertainty in the platescale and field distortion. As another cross-check to the accuracy, we applied this procedure to an additional epoch of F220M data from 1998 February 11 (which had a significantly different roll angle), and the two epochs agreed well within assigned uncertainty, indicating that our 35 mas estimate for error in the diffraction-spike centroid method is indeed quite conservative. Complete results appear in tabulated and graphical forms in Table 2 and Figure 3. HST imagery is rarely used to perform relative astrometry at this level of precision, and it is possible that unknown calibration errors may corrupt these results. Since the registration of the infrared and visible images is critical for interpreting the emission morphology, we are hopeful that independent confirmation of these results can be obtained.

#### 2.3.3 Proper motion of CIT 6

Because the WFPC2 and NICMOS images were taken two years apart, it is important to consider the possibility that the relative astrometry is corrupted due to motion of CIT 6 with respect to the reference star. Bastian et al. (1988) documented the binary separation of these sources as 10.4″$`\pm `$0.4″(PA 18°$`\pm `$1°) at epoch 1955.3 using Palomar Sky Survey (red) plates. The red WFPC2 data from 1996.3 shows a separation of 11.14″$`\pm `$0.02″(PA -20.3°$`\pm `$0.1°). According to these results, the proper motion of CIT 6 with respect to the reference star is 0.021″$`\pm `$0.014″yr-1. This implies a 0.040″$`\pm `$0.027″shift between the WFPC2 and NICMOS observations. This is at the same level as the error on the K band centroid measurement ($`\pm `$0.035″). Importantly, if one were to apply this correction it would cause the infrared centroids in Figure 3 to shift away from the visible centroid towards the North-East.

## 3 Discussion

Because the northern component of the Keck 2.2µm images (see Figure 1) dominates the flux, the centroid of the NICMOS images should lie close to this location. Based on the HST relative astrometry, it appears that the northern component in the IR images can be associated with the northern component in the red WFPC2 images, although they do not appear to be exactly co-spatial. Interestingly, this places the faint, southern IR component approximately halfway between the two visible components. Note that the small proper motion discussed in §2.3.3 would cause the IR centroids to appear further to the North-East, away from the average position of the visible components, strengthening the above interpretation. One interpretation of the two visible components is that they arise from a bipolar reflection nebula with low density poles allowing scattered light to escape the dust shell. This is supported by the linear polarization of the red emission, with a PA of $``$75°(Trammell et al. 1994). The dominance of the northern component in the IR indicates that the northern pole is inclined towards us. In this scenario, the faint, southern IR component could be thermal emission from the central “disk” of material which is responsible for channeling the stellar flux towards the bipolar lobes.

Unfortunately, this appealing model does not account for the blue WFPC2 observations. In fact, the southern component is brighter in all the WFPC2 images. However, the source of blue light is likely to be different from the source of red emission (Alksnis et al. 1988). This contention is most securely based on multi-wavelength photometry which shows that the red and infrared flux vary at the $``$626 day period of the carbon star, while the blue flux level is roughly constant (Alksnis 1995). In addition, the blue light has a polarization angle ($``$0°) roughly orthogonal to that seen in the red and infrared (Cohen & Schmidt 1982; Trammell et al. 1994).

One scenario able to reconcile these disparate observations is to propose that the blue lobe is not being illuminated the carbon star, but by a blue companion which lies in the vicinity of the southern lobe. Although it may appear somewhat ad-hoc, this model is able to qualitatively explain the odd cometary shape of the bluest WFPC2 image and the different blue-red polarization angles. The apparent location of the southern (blue) HST component would be largely defined by the mass-loss profile of the carbon star, while the angle of linear polarization in the blue would be set by the separation vector between the companion and the scattering lobe. Note that the polarization direction of the Hα emission observed by Trammell et al. (1994) matches the blue polarization, not the red. This is consistent with the idea that the blue companion is producing the line emission which is scattering into our line of sight along the same trajectories as the blue continuum.

### 3.1 A Dusty Disk?

Generally, scattered photospheric light should have a higher color temperature than thermal emission from dust. Hence, the above model can be crudely tested by calculating the color temperature of the northern and southern IR components. Towards this goal, all the epochs of Keck data have been averaged together, where image registration was estimated by maximizing the cross-correlation; these coadded images appear in Figure 4. It is clear that the southern component has relatively more flux (compared to the northern component) at 3.1µm than at 2.2µm, consistent with the southern component being redder. These images have been projected onto an axis at PA 30°, and the summed profiles appear in Figure 5. Note that the larger FWHM of the 3.1µm peak is due to the lower spatial resolution, and the profiles have been normalized by the total flux under each curve rather than the peak intensity.

To minimize bias due to the wavelength-dependent spatial resolution, two photometric apertures have been defined to include most of the flux in the northern and southern IR components. These apertures appear in Figure 4 and their projected boundaries along PA 30°are marked on Figure 5. At 2.2µm (3.1µm), the southern aperture contains $``$8% ($``$18%) of the total flux in the two apertures. Taranova & Shenavrin (1999) report K and L band magnitudes for CIT 6 during the period of these Keck observations, with a mean K mag ($`\lambda _0=2.2\mathrm{\mu m}`$) of 1.83, and L mag ($`\lambda _0=3.5\mathrm{\mu m}`$) of -1.10. We linearly interpolate (in wavelength–magnitude space) to arrive at an approximate 3.1µm magnitude of 0.22. Using these average fluxes we can now estimate the average color temperature of the northern and southern components separately.

The color temperatures based on the 2.2µm and 3.1µm images are 870 K and 600 K for the northern and southern components respectively. No color corrections were applied when converting magnitudes to physical flux units, but any corrections will modify the component temperature estimates in the same way. As a comparison, the mean color temperature based on J and K-band magnitudes is $``$1000 K, according to Taranova & Shenavrin (1999). Although the computation above is clearly an approximation, the southern IR component has been found to be demonstrably redder than the northern one.

However, the lower color temperature of the southern component does not guarantee that this is thermal emission from dust. The relative contributions of reddened scattered light and thermal emission are difficult to untangle in optically thick envelopes (e.g., Weigelt et al. 1998; Monnier et al. 1999; Tuthill et al. 2000a). The fact that the southern IR component is located roughly between the visible lobes and has the lowest color temperature, is certainly consistent with it being a dusty “disk” of material surrounding the carbon star.

### 3.2 Time-variability

The dominant (northern) IR component clearly changes shape on the time scale of years, as can be seen from examination of Figure 1. Most notably, the northern IR component has evolved from a nearly symmetric elongation in 1997 to a very asymmetric structure extending West-Northwest. These changes are observed at both wavelengths and show a consistent evolution in time. Under the above hypothesis that this component represents scattered light or thermal emission for warmed dust, this time evolution can be understood if dust formation closer to the star is changing the illumination pattern. These observations represent the first direct detection of evolution of the inner dust shell of CIT 6. The inner dust formation zone around another nearby carbon star, IRC+10216, has also been observed to be quite chaotic and dynamic, evolving on similar time scales (Tuthill et al. 2000b). However in IRC +10216, discrete dust blobs are observed to move away from central souce, while only slight shifts in the emission pattern are seen for CIT 6. This difference can be understood if the central star of CIT 6 is more obscured than for IRC +10216, and also by considering that CIT 6 is thought to be be $``$4 times more distant than IRC +10216.

If illuminated by scattered light, this changing morphology of the northern IR component could explain the observed changes in polarization angle of the red and infrared light. Since 1970, the polarization angle has varied by $`\pm `$20°, perhaps related to the opening angle of the northern lobe, as can be seen by the polarization vectors placed on Figure 2. As the center of scattered light changes in this lobe, the corresponding polarization angle would as well. Coordinated spectropolarimetry and high spatial resolution imaging could test this hypothesis further.

### 3.3 Alternate models

Unfortunately, these new data do not preclude other interpretations. Alternate possible dust shell geometries include the northern component being the star itself (or a binary), with the southern component arising from scattered radiation. As for the case of carbon star IRC +10216, it is not obvious if the star itself is visible in existing near-infrared images, or if it remains obscured (see discussions in Haniff & Buscher 1998; Weigelt et al. 1998; Tuthill et al. 2000a). Yet another possibility is that the elongated and time-variable morphology of the northern IR component could be due to a close binary.

A direct measure of the polarization percentage of both IR components would more definitively determine which (if either) are due to scattered radiation, although such observations are technically difficult. Alternate theories for the blue excess, either due to a mistake in the spectral type identification (Trammell et al. 1994) or circumstellar shocks related to mass-loss (Cohen & Schmidt 1982), do not seem able to explain either the odd cometary appearance or the polarization properties without positing an additional bipolar nebula at right-angles to the first. Given the IR morphology, the wide binary hypothesis seems the best hypothesis of those presently considered for explaining the peculiar properties of the blue light, but stronger evidence is needed to confirm it.

## 4 Conclusions

The dust shell of carbon star CIT 6 has been examined with unprecedented spatial resolution ($`\stackrel{<}{_{}}`$50 mas) using aperture masking at the Keck-I telescope. Serendipitously, a nearby ($``$11″) companion allowed high precision relative astrometry from recent HST WFPC2 and NICMOS images. Multi-wavelength analysis suggests that CIT 6 is surrounded by a bipolar dust envelope with a complex morphology, consisting of at least two components. The dominant compact feature in the near-infrared is tentatively identified as the northern scattering lobe of the reflection nebula, although the possibility that it is the carbon star itself can not be ruled out. In addition, a fainter IR component with very red colors has been detected between the two visible light components, perhaps due to thermal emission from the central “disk” of the bipolar dust envelope. Excess blue emission in the southern optical lobe is ascribed to an obscured hot companion, responsible for the strong wavelength-dependence of the (visible) polarization properties and “cometary” scattering feature in the blue WFPC2 image. Under this model, changes in the red and infrared polarization angle should be linked to observed variations in the near-IR emission morphology.

Other recent high resolution observations of evolved stars have also revealed strong deviations form spherically symmetry. Red supergiant VY CMa (Monnier et al. 1999) and carbon star IRC +10216 (Haniff & Buscher 1998; Weigelt et al. 1998; Tuthill et al. 2000a) also show evidence for bipolar nebulae on large scales (as seen in scattered light) and complicated dust shells in the innermost regions (from IR observations). The dominant physical mechanism for breaking spherical symmetry and producing the bipolar mass-loss profile remains unidentified.

The authors would like to recognize Dean Hines, Gary Schmidt, and Susan Trammell for sharing results of their independent studies in advance of publication and for stimulating conversations regarding the CIT 6 dust shell. In addition, we thank Eddie Bergeron and Glenn Schneider for advice regarding relative astrometry with HST. The image reconstructions presented here were produced by the maximum-entropy mapping program “VLBMEM,” written by Devinder Sivia. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and NASA’s Astrophysics Data System Abstract Service. The data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. These results are based also on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. This work is a part of a long-standing interferometry program at U.C. Berkeley, supported by the National Science Foundation (Grants AST-9321289 and AST-9731625) and NASA. JDM was partially supported by a Postdoctoral Fellowship from the Harvard-Smithsonian Center for Astrophysics.