text:
# Optical spectra of the carbon-oxygen accretion discs in the ultra-compact X-ray binaries 4U 0614+09, 4U 1543-624 and 2S 0918-549

## 1 Introduction

Low-mass X-ray binaries are systems in which a neutron star or black hole accretes from a low-mass companion. Most systems have orbital periods of hours to days and are consistent with the scenario (van den Heuvel, 1983) in which the donors are main sequence or evolved, hydrogen-rich, stars. A few ultra-compact systems have orbital periods below an hour and are so compact that the donor stars cannot be main sequence stars, but instead must be hydrogen poor (e.g. Verbunt & van den Heuvel, 1995).

From X-ray spectra of five X-ray binaries, including the persistent sources 4U 0614+09, 4U 1543-624 and 2S 0918-549 Juett et al. (2001) inferred an enhanced neon/oxygen ratio, which they interpreted as being local to the systems. The similarities between these three systems and the other two, which are known ultra-compact X-ray binaries, led Juett et al. (2001) to conclude that these systems all have ultra-short orbital periods and to propose that their donor stars originally were carbon-oxygen or oxygen-neon-magnesium white dwarfs. Recently, Juett & Chakrabarty (2003) reported further X-ray spectroscopy for 4U 1543-624 and 2S 0918-549, confirming their earlier findings. Yungelson, Nelemans & van den Heuvel (2002) argued, based on mass-transfer stability arguments, that the donor stars in ultra-compact X-ray binaries that have formed from white dwarf – neutron star binaries should be low-mass white dwarfs ($`M_{\mathrm{donor}}0.45\mathrm{M}_{}\text{ }`$). Combining binary evolution constraints and white dwarf interior studies, Yungelson et al. (2002) concluded that these systems could be brought into one unifying scheme in which the systems were descendants of binaries consisting of a so called hybrid white dwarf (carbon-oxygen (CO) core with thick helium mantle, Iben & Tutukov, 1985) and a neutron star. The systems came into contact by angular momentum loss due to gravitational wave emission, and quickly evolved from periods of a few minutes to typical periods of tens of minutes. After cooling for several Gyr the interior of a CO white dwarf crystallises and, due to differential gravitational settling, chemical fractionation will occur (e.g. Hernanz et al., 1994). At periods above 10 minutes neon enriched layers could be exposed (see Yungelson et al., 2002). At these orbital periods the only alternative donor stars are degenerate helium stars or hydrogen poor remnants of stars that started mass transfer at the very end of the main sequence (for the latter, see Podsiadlowski et al., 2003), which would consist mainly of helium.

To further test the possible ultra-compact nature of 4U 0614+09, 2S 0918-549 and 4U 1543-624, we obtained optical spectra for these sources, because ultra-compact systems are expected to be hydrogen deficient, and even more, neon-rich donors stars are also expected to be helium deficient. Previous optical spectra of 4U 0614+09 (Davidsen et al., 1974; Machin et al., 1990) indeed show no signs of the classical accretion disc hydrogen emission lines but have a relatively low S/N ratio.

## 2 Observations and reduction

Spectra were taken with the FORS2 spectrograph on UT4 of the 8m Very Large Telescope on Paranal in Chile. For each object we took spectra both with the 1400V and 600RI holographic grisms, with a 1” slit, using 2x2 on-chip binning. This setup resulted in coverage of 4620 – 5930 Å with mean dispersion of 0.64 Å/pix for the 1400V spectra and 5290 – 8620 Å with mean dispersion of 1.63 Å/pix for the 600RI spectra. A log of the observations is given in Table 1.

Data reduction was done using standard IRAF111IRAF is distributed by the National Optical Astronomy Observatories tasks. The bias was removed using the overscan region of the CCD, after which the images were flatfield corrected using the standard calibration plan flatfields. Spectra were extracted using optimal extraction (Horne, 1986) with the apall task. Arc lamp spectra were extracted from the same place on the CCD. The 1400V wavelength calibration was obtained using the positions of 17 lines, giving a root-mean-square scatter of 0.05 Å in fitting a fourth-order Lagrangian polynomial. The 600RI wavelength calibration uses 40 lines and gives a root-mean-square scatter of 0.15 Å for a fourth-order Lagrangian polynomial.

The spectra were flux calibrated, using the nearest two flux standard stars that were available from the VLT archive. Since these were taken many days from our observations the flux calibration only provides a very rough absolute flux calibration, but it does give a reasonable estimate of the continuum shape. The reduced spectra were subsequently imported in the MOLLY package for further analysis. For 4U 0614+09, where the S/N ratio is the highest, we removed the telluric absorption features by dividing the flux calibrated spectrum by a template of the absorption. The template was constructed from the spectrum of a bright star that was also in the slit as follows. We fitted a third order cubic spline to line-free regions of the continuum of this spectrum and used the fit to normalise the continuum to unity. The value of all pixels in the spectrum was set to 1, except for the wavelength rages 6865–7700 and 8085 – 8265 Å, where strong telluric features are present.

For each object all spectra were combined and averaged to obtain one final spectrum.

## 3 Analysis and interpretation

### 3.1 Identification of the lines

In Fig. 1 we show the combined spectra. In the figure we also show the positions of the most common hydrogen emission lines seen in classical X-ray binaries as solid vertical lines (H$`\alpha `$ at 6563 Å and H$`\beta `$ at 4861 Å). We also show the positions of several helium emission lines observed in optical spectra of some interacting binaries (e.g. GP Com, see Marsh et al., 1991). The dashed lines show He i lines at 4713, 4921, 5015, 5876, 6678, 7065 and 7286 Å. The dotted line shows the position of the He ii 4686 Å line often seen together with hydrogen emission. From Fig. 1 it is clear that the objects do not show strong lines of hydrogen or helium. Upper limits on the equivalent widths (EWs) of these lines are shown in Table 2 (where we give the value of $``$EW to avoid minus signs). There is a hint of absorption around H$`\beta `$ in 4U 1543-624 and possibly emission around He i 4713 in 2S 0918-54.

The 4U 0614+09 spectrum shows clear, though sometimes weak, emission lines. We propose that these are due to partially ionised carbon and oxygen, which have many lines in the optical region. In order to identify the lines we calculated a simple emission profile for a 30-70 mixture (by number) of carbon and oxygen at a temperature of 27,000 K, using an LTE emission line model, as described in Marsh et al. (1991). We assume a particle density of 3 $`\times `$ 1013 cm-3 and a line of sight through the medium of 107 cm. The only update is that we consider ions up to C v and O v and use an atomic line list compiled from the National Institute of Standards and Technology222Version 2.0, http://physics.nist.gov/cgi-bin/AtData/lines\_form and the Atomic Line List333http://www.pa.uky.edu/$``$peter/atomic/. The model is smoothed with a 12 Å Gaussian kernel, which at these wavelengths corresponds to a $``$600 km s-1 velocity width, *if* it is associated with kinematics. This model is not a physical model of the accretion disc, as photoionisation is surely important. However, we believe it correctly identifies the lines in the observed spectrum. Almost all features are blends of multiple lines, most of them from C ii, C iii and O ii, while the strong feature at 5585 is due to O iii. Around 5810 C iv might be present. In Table 3 we give the positions of the strongest features and the lines that may contribute to them. We also list the measured equivalent widths of the different features (again negative). The entries within parentheses are either weak lines or uncertain identifications. The 7240 and 8220 features could be affected by the telluric correction.

The lower two spectra in Fig. 1 are of 4U 1543-624 and 2S 0918-549, rebinned to 4 Å to enhance the S/N ratio. In Table 3 we also list the equivalent widths of the features in these spectra. The spectrum of 4U 1543-624 looks like a scaled down version of the 4U 0614+09 spectrum, except maybe for the strength of the C iii and O ii complexes around 4650 and 4700 Å. All features, except the carbon features at 5280, 6070 and 6150 Å are detected. For 2S 0918-549 the case is less clear, mainly because the lower signal. The features around 4650, 4700 and 6730 Å are present, but most of the other features are undetected.

### 3.2 Spectral shape and reddening

In Table 2 we also list the measured spectral slope of the three objects, assuming a power law ($`F_\lambda \lambda ^\alpha `$), and the temperatures of the best fitting black body spectra ($`T_{\mathrm{BB}}`$). The spectrum of 4U 1543-624 is significant bluer than the other two. For 4U 0614+09, the continuum is well fit with $`\alpha =1.67`$ for $`\lambda 5100`$ Å, but is significant flatter below. This flattening is not seen in the spectrum of Machin et al. (1990), and thus might be due to the flux calibration.

We also list several estimates or limits on the interstellar reddening of the objects. We have two possible reddening indicators available in the spectra: the Na D lines and the diffuse interstellar band (DIB) at 5780 Å. In Table 2 we list the measure equivalent widths of these features and the implied reddening according to Munari & Zwitter (1997) and Herbig (1993) respectively. The Na D lines only give rough lower limits, as the EWs level off above $`E(BV)`$ 0.5 due to saturation. The high EWs measured here would then be the result of multiple components in the lines, as can indeed be seen in the 4U 0614+09 spectrum. In the Table we also list the maximum reddening according to the Schlegel et al. (1998) dust maps and the implied reddening for the hydrogen column found in the spectral modelling of the X-ray spectra (Juett et al., 2001). These last values are upper limits, as the local O and Ne absorption implied by the spectra would artificially enhance the hydrogen column in the fits to the X-ray spectra (e.g. Juett et al., 2001). Although the above gives a barely consistent picture, we conclude that the reddening of these objects is probably close to the maximum in the Galaxy in their directions.

### 3.3 Origin of the line emission

Although detailed modelling of the observed spectra is beyond the scope of this paper, we briefly discuss the possible origin of the line emission. The spectral shapes of the continua, corrected for the reddening (using $`E(BV)`$ = 0.6, 0.4 and 0.6 for 4U 0614+09, 4U 1543-624 and 2S 0918-549 respectively) are consistent with a black body of temperature $``$20,000 K. The fact that these objects are persistent X-ray sources, could point at origin of the lines in either the irradiated disc or the irradiated donor star, although in that case more than doubly ionised species might be expected to be present. Only phase resolved spectroscopy will conclusively tell, but the observed width of the lines ($`600`$ km s-1) might suggest the inner parts of the accretion disc.

## 4 Discussion

Previous spectra have been interpreted in a different way to that presented here, but were hampered by quite low S/N ratio. For instance Machin et al. (1990), interpreted the strongest lines near 4650 Å and 5590 Å, as the Bowen blend and possibly O i at 5577 Å. They already commented on the intriguing absence of the usual He II 4686 Å line, which casts doubt on the interpretation of the emission at 4650 Å as the Bowen blend, since the Bowen mechanism is driven by helium. The total absence of any sign of helium at any of the positions indicated in the plots strongly argues against (much) helium in the system.

As discussed in the Introduction, Juett et al. (2001) suggested the donor stars originally were carbon-oxygen or oxygen-neon-magnesium white dwarf, while Yungelson et al. (2002) argued they were hybrid white dwarfs. Our spectra thus give further evidence for the interpretation of 4U 0614+09, 2S 0918-549 and 4U 1543-624 as ultra-compact X-ray binaries with carbon-oxygen white dwarf donors. One very interesting prospect is the determination of the carbon to oxygen ratio in the transferred material from detailed modelling of the disc spectrum, which would give an unprecedented view into the interior of a white dwarf and possibly could be used to constrain the rate of the $`{}_{}{}^{12}C(\alpha ,\gamma )^{16}O`$ reaction (see for a discussion and references Straniero et al., 2003).

Further evidence for the ultra-compact nature of 4U 0614+09 and 2S 0918-549 comes from their absolute visual magnitude. Using the magnitudes and reddening of the objects from the low-mass X-ray binary catalogue (Liu et al., 2001) and the distances of these systems of $`<`$3 and 4.2 kpc respectively from type I X-ray bursts (Brandt et al., 1992; Cornelisse et al., 2002) we find absolute magnitudes of $`>`$5.4 and 6.9 for 4U 0614+09 and 2S 0918-549. These are very faint for X-ray binaries and, from the van Paradijs & McClintock (1994) relation between absolute magnitude, X-ray luminosity and orbital period, suggest periods well below one hour. The derived X-ray luminosities of $`0.01L_{\mathrm{Edd}}`$ for 4U 0614+09 (Ford et al., 1996) and $`0.003L_{\mathrm{Edd}}`$ for 2S 0918-549 (Jonker et al., 2001), assuming these are indicative of the average mass accretion rate, give mass transfer rates of $``$3 $`\times `$ 10-10 and $``$9 $`\times `$ 10-11 $`\mathrm{M}_{}`$ yr-1 for 4U 0614+09 and 2S 0918-549, assuming accretion of helium in the type I X-ray bursts (see below). Comparing these rates with Fig. 15 of Deloye & Bildsten (2003), suggests orbital periods of 15-20 and 20-30 minutes for 4U 0614+09 and 2S 0918-549.

Finally we note that our findings pose an interesting question concerning type I bursts. As discussed by Juett & Chakrabarty (2003), the bursts observed in 4U 0614+09 (Swank et al., 1978; Brandt et al., 1992) and 2S 0918-549 (Jonker et al., 2001) are all short bursts, which are believed to be caused by hydrogen and/or helium. They suggest that possibly the donor stars still have non-negligible hydrogen fraction, which would have evolved from binaries that start mass transfer close to the end of the main sequence (Nelson et al., 1986; Podsiadlowski et al., 2002). However, these are expected to consist mainly of helium, rather than carbon and oxygen. The alternative Juett & Chakrabarty (2003) suggest, which seems the only remaining option in light of the lack of helium lines in the optical spectrum, is that spallation of the carbon and oxygen nuclei at the impact onto the neutron star (cf. Bildsten et al., 1992) turns them into helium (or possibly hydrogen), which subsequently triggers the burst.

## 5 Conclusions

We presented optical spectra of the three suspected ultra-compact X-ray binaries 4U 0614+09, 4U 1543-624 and 2S 0918-549. The spectra show no sign of hydrogen or helium emission lines. We identify the observed features as lines of C ii, C iii, O ii and O iii. This is in agreement with the interpretation of these sources as ultra-compact X-ray binaries and the expectation that the donor stars in these objects originally were hybrid white dwarfs that have lost most of their mass and now consists mainly of oxygen and carbon. These spectra are thus of accretion discs almost purely made out of oxygen and carbon.

## Acknowledgments

We thank the referee Janet Drew for comments that improved the paper. We are thankful to Peter van Hoof and the National Institute of Standards and Technology for compiling the atomic line lists we use. We further thank Lars Bildsten for stimulating discussions. GN acknowledge the hospitality of the Kavli Institute for Theoretical Physics. This work was supported by the National Science Foundation under grant PHY99-07949.


format: text

text:
# The rotational broadening and the mass of the donor star of GRS 1915+105

## 1 Introduction

The microquasar GRS 1915+105 remains one of the most exotic and variable objects in the X-ray, IR and radio spectrum bands since its discovery (Castro-Tirado et al. 1994). A major accomplishment has been the implied connection between the jet and the accretion disk as interpreted through modelling of its X-ray variability (Belloni et al. 1997; Klein-Wolt et al. 2002). Thereafter, it is evident that crucial insight into the accretion physics of jet formation can be derived through intense theoretical and observational work. Fundamental for the modelling but also for the evolutionary history of GRS 1915+105 are the binary system parameters.

The binary model of the microquasar GRS 1915+105 was elusive, until recently, despite many observational efforts to search for the donor star (Castro-Tirado et al. 1996; Eikenberry et al. 1998; Mirabel et al. 1997; Martí et al. 2000; Harlaftis et al. 2001; Greiner et al. 2001a). Greiner et al. (2001b) were able to obtain phase-resolved IR spectra using the VLT-Antu with the ISAAC IR spectrograph and derive the orbital period of the binary system from the detected Doppler-shifted 12CO and 13CO absorption bands. The deduction of the masses of the binary, using the mass function equation,

$$f(M_\mathrm{x})=\frac{PK_\mathrm{d}^3}{2\pi G}=\frac{M_\mathrm{x}\mathrm{sin}^3i}{(1+q)^2},$$

requires the orbital period $`P`$, the semi-amplitude of the radial velocity curve of the donor star around the black-hole $`K_\mathrm{d}`$, the binary inclination $`i`$ and the mass ratio $`q=M_\mathrm{d}/M_\mathrm{x}`$, where $`M_\mathrm{d}`$ is the mass of the donor star and $`M_\mathrm{x}`$ the mass of the black-hole. Indeed, the orbital period and the semi-amplitude of the radial velocity curve of the donor star around the black-hole have been determined by Greiner et al. (2001b; $`P=33.5\pm 1.5`$ days and $`Kw_\mathrm{d}=140\pm 15`$ km s-1). The remaining unknown parameters are the binary inclination and the mass ratio. In fact, it is generally the binary inclination which dominates the uncertainty in the mass determination of the black hole via $`f(M_\mathrm{x})\mathrm{sin}^3i`$. The general method in deriving the inclination of X-ray binaries and cataclysmic variables is to model the ellipsoidal modulations of the companion star at infrared wavelengths. Luckily, the orbital inclination of GRS 1915+105 is deduced with unprecedented accuracy from the orientation of the jets ($`i=70^{}\pm 2^{}`$; Mirabel & Rodriguez 1994; $`i=66^{}\pm 2^{}`$ by Fender et al. 1999). The final parameter that is still undetermined for GRS 1915+105 is its mass ratio. The mass ratio of the binary system is related to the rotational broadening of the donor star using the formula

$$\frac{\upsilon \mathrm{sin}i}{K_\mathrm{d}}=0.46\left[(1+q)^2q\right]^{1/3}$$

where $`\upsilon \mathrm{sin}i`$ is the rotational broadening in the line-of-sight, assuming Roche lobe mass-transfer, spherical geometry for the star and tidal locking to the primary star (Gies and Bolton 1986; Horne, Wade, Szkody 1986). Here, we further analysed the same IR spectra from VLT (Greiner et al. 2001b) in order to deduce the rotational broadening of the donor star from the width of the photospheric absorption lines (mainly 12CO, 13CO; see Greiner et al. 2001a).

## 2 Observations, numerical technique and deduction of mass ratio

The GRS 1915+105 spectra, together with a KIII spectrum (HD 138185), were obtained between April and August 2000 with the ISAAC spectrograph on the Antu-VLT covering the range 2.2896–2.4131 $`\mu `$m, at a pixel resolution of 15.4 km s-1 (Greiner et al. 2001b). In order to measure the width of the photospheric absorption lines, we apply a numerical technique according to which we subtract from the Doppler-corrected spectrum of GRS 1915+105 a KIII template spectrum. In detail

* we produce the average spectrum after the individual spectra are Doppler-corrected by applying the radial velocity solution (Greiner et al. 2001b)
* we apply to the KIII template spectrum a rotational profile with a limb darkening coefficient of 0.5 (Gray 1992) and a veiling factor $`f`$ which scales the depth of the absorption lines in order to simulate the GRS 1915+105 absorption lines
* the difference spectrum represents mainly the accretion disk spectrum with all the photospheric absorption bands arising from the KIII donor star of GRS 1915+105 having been subtracted
* we perform a $`\chi ^2`$-test between the difference spectrum and its smoothed version (by applying a Gaussian profile of a large FWHM=450 km s-1 to make it flat) which is a measure of how good the absorption lines of the donor star have been subtracted for each rotational profile and scaling factor $`0<f<1`$ applied
* we iterate by changing the rotational profile width and the depth of the absorption lines - using a scaling factor $`f`$ representing the the donor star line flux - until we reach $`\chi ^2`$ values moving along a parabolic surface.

The minimum $`\chi ^2`$ on the parabola ($`\upsilon \mathrm{sin}i,\chi ^2`$) gives the corresponding width of the rotational profile. We deduce a $`\upsilon \mathrm{sin}i=26\pm 3`$ km s-1 and $`f=14\pm 3`$ %, at the 90% significance level, where $`f`$ is the light contribution from the donor star as estimated from its 12CO/13CO absorption bands. As a validation test, we performed the same iterations but now for the individual spectra rather than the average Doppler-corrected spectrum and took the average $`\upsilon \mathrm{sin}i`$ \- equal to $`26\pm 5`$ km s-1 \- from the individual estimates. The mass ratio deduced is, only weakly sensitive, on the spectral type and veiling factor $`f`$. For more detailed discussion of the technique see Harlaftis et al. (1999 and references therein). Fig. 1 shows, from bottom to top, the KIII template (with lines multiplied by the veiling factor of 0.14), the average Doppler-corrected spectrum of GRS 1915+105 (16 spectra), the difference spectrum (thus, representing the disk spectrum) and the difference spectrum between the template and the average Doppler-corrected spectrum of GRS 1915+105 at binary phase $``$0.75 (same average as in Greiner et al. 2001b). The KIII spectrum has a ripple in the continuum arising from variations in the star brightness between integrations (i.e. at different detector positions). There is a residual blue-shifted CO-band or other unidentified emission in the disk spectrum from phases other than 0.75 as well as two unidentified absorption features around 2.2962 and 2.2988 $`\mu `$m.

## 3 Discussion

We discuss now the various assumptions involved in the derivation of the mass ratio from the rotational broadening. The absorption vibration-rotation bands 12CO and 13CO consist of $`4500`$ lines at 3500 K (Kunde 1977). However the total damping width is of the order of $`1\times 10^4`$ Å , the thermal width does not exceed 0.04 Å at 4000 K, and the macroturbulence width is $`5`$ km s-1 in giant stars (Gray 1992). Other possible broadening factors other than rotation is pressure-broadening which is not significant in giant stars. Thus, the fact that there is a measurable velocity dispersion in the line profiles indicates that it is dominated by the rotational broadening profile. Shahbaz (2003) has developed a new, more accurate method involving synthetic spectra from the Roche lobe (thus avoiding assumptions about a spherical symmetric star with a limb darkening coefficient). However, simulations show that the standard method compared to his more accurate treatment is still compatible within 3$`\sigma `$. The assumption of a spherical star in the ($`\upsilon \mathrm{sin}i,q`$) relation has been shown to introduce an error in the evaluation of mass ratio $`q`$ of $`<5`$ % (Marsh et al. 1994) which is much less than the uncertainty of our measurement ($`50`$ %), hence with a negligible effect. The assumption of synchronous rotation ($`P_{\mathrm{orbit}}=P_{\mathrm{rot}}`$) has been found to be valid for binaries as wide as 80 days even though giant stars may show a weaker synchronism (10% difference from the orbital period) than late-type stars (Giuricin et al. 1984). The criterion for tidal locking is $`R_\mathrm{d}/\alpha >0.125`$, where $`R_\mathrm{d}`$ is the radius of the donor star and $`\alpha `$ the binary separation, assuming dynamical tide alone (Zahn 1977; Zahn 1992) which indeed holds for GRS 1915+105 ($`R_\mathrm{d}/\alpha =0.176`$). Involving other dissipative processes such as hydrodynamical meridional flow on the star lowers the limit for the fractional star radius, thus yielding much longer periods (Tassoul 1987; Tassoul & Tassoul 1992; see also for tidal interaction models Witte & Savonije 2002, Claret & Cunha 1997, Verbunt & Phinney 1995). The synchronisation timescale $`t(R_\mathrm{d}/\alpha )^6q^2=4.6`$ million years (Giuricin et al. 1984) which is well within the nuclear evolution timescale of the donor star. Synchronisation is more rapid than circularisation of the orbit (Hall 1986; see also Lecar et al. 1976), and a $`\chi _\nu ^2`$-test on the radial velocity data (Greiner et al. 2001b) does not justify an eccentric orbit ($`e<0.006`$). Therefore, we conclude that the orbit is circular.

The constraint on the eccentricity of the radial velocity fit also indicates that there is no line irradiation affecting the measurement of the $`K`$ semi-amplitude of the radial velocity of GRS 1915+105. For example, a 10% correction is applied in the value of the $`K`$ semi-amplitude (reduced) from the NaI doublet absorption lines of the donor star in the eclipsing dwarf nova IP Peg from a measured eccentricity of the radial velocity fit of ($`e=0.089`$; Martin et al. 1989). There has been no evidence for irradiation affecting the photospheric absorption lines in any quiescent black-hole X-ray transients with extreme mass ratios (e.g. Harlaftis et al. 1996). However, outbursts can irradiate the donor star and affect significantly the radial velocity amplitude, as for example with the FIII-IV sub-giant donor star in GRO J1655-40 (Shahbaz et al. 2000) which has a binary period of 2.6 days and a large mass ratio of 0.42 (Shahbaz 2003). The spectra of GRS 1915+105 were obtained during relatively quiet states (no outburst or flare) and in addition its binary orbit is long with a small cross section of the donor star for any irradiation to affect it. In fact, using the relation (2) from Phillips et al. (1999) we find the effect to be only 2.5% compared to the 11% uncertainty in $`K_\mathrm{d}`$ for GRS 1915+105.

The mass ratio thus derived from the rotational broadening width is $`q=0.058\pm 0.033`$ (see second relation in the Introduction). The masses can now be determined from $`(P,K,i,q)=`$($`33.5\pm 1.5`$ d, $`140\pm 15`$ km s-1, $`66^{}\pm 2^{}`$, $`0.058\pm 0.033`$) giving $`M_\mathrm{x}=14.0\pm 4.4M_{}`$ and $`M_\mathrm{d}=0.81\pm 0.53M_{}`$ (see first relation in the Introduction). The large black-hole mass uncertainty is dominated by the uncertainty in the semi-amplitude radial velocity of 15 km s-1. The contributions from the period uncertainty and mass ratio are about equal but only a few % whereas that of the inclination uncertainty is practically zero. The large mass ratio uncertainty of 50% is due to the large uncertainty in both the rotational broadening and the semi-amplitude of the radial velocity ($`10`$% for each parameter). Having made the above points for the assumptions involved we now discuss the ambiguity that may be reside in the inclination value. The jet axis may not be perpendicular to the binary plane though a constraint is placed to be less than $`25\mathrm{°}`$ in GRS 1915+105 (Maccarone 2002). Indeed, in one of the microquasars, namely GRO J1655-40, the binary inclination is $`70\mathrm{°}`$ (Greene et al. 2001) whereas the jet inclination is found at $`85^{}`$ (Hjellming & Rupen 1995). The proper motions of the plasmoids ejected from GRS 1915+105 have been found to vary in apparent velocity and this was attributed to a change in the jet velocity (Rodriguez & Mirabel 1999; Fender et al. 1999). However, it must be noted that the alternative explanation is that the jet axis may be precessing. Indeed, Rodriguez & Mirabel (1999) mention a change by $`10^{}`$ of the ejection axis within one month as well as a 10–20 % difference in the proper motion measured by Fender et al (1999). However, Maccarone (2002) concludes, on the basis of close agreement between the dynamical estimate and the resonance interpretation of the observed X-ray QPO frequencies for the black hole mass, that a large misalignment between the jet axis and the disk axis is ruled out in GRS 1915+105.

The above discussion must be taken under consideration when using the adapted inclination of $`i=66^{}\pm 2^{}`$ (Fender et al. 1999) for a distance of $`11.2\pm 0.8`$ kpc. Assuming a $`10^{}`$ uncertainty in the inclination due to possible jet axis precession, the masses become $`M_\mathrm{x}=11.6\pm 3.3M_{}`$ and $`M_\mathrm{d}=0.68\pm 0.43M_{}`$ for an inclination of $`i=76^{}`$ ($`M_\mathrm{x}=16.9\pm 5.9M_{}`$ and $`M_\mathrm{d}=0.98\pm 0.65M_{}`$ for $`i=56^{}`$). However, extensive radio monitoring of GRS 1915+105 over the last 10 years indicates no measurable constant precession, thus giving ground to the interpretation of varying jet ejection rather than a jet axis precession. For direct comparison on the improvement of the mass estimate using the mass ratio determined here, we deduce a $`M_\mathrm{x}=12.9\pm 4.04M_{}`$ compared to $`M_\mathrm{x}=14.0\pm 4.0M_{}`$ for $`i=70^{}\pm 2^{}`$ as adapted by Greiner et al. (2001a).

Single KIII stars have negligible $`\upsilon \mathrm{sin}i`$. Indeed, a K1III and a K5III giant star have a mass of 2.3, 2.2 $`M_{}`$, a radius of 11, 28 $`R_{}`$ and a $`<\upsilon \mathrm{sin}i>`$ of 2.5, and less than 1.5 km s-1, respectively (Gray 1992). Greiner et al. (2001a) identified the luminosity class of the donor star in GRS 1915+105 as a giant star mainly based on the equivalent width ratio of the 12CO/13CO lines. Using the Eggleton formula (1983) to determine the Roche lobe size and hereafter the density of the lobe-filling star, we derive $`R_\mathrm{d}=19R_{}`$ and $`\rho =2\times 10^4`$ g cm-3, using an orbital period of $`P=33.5`$ days and a mass ratio of $`q=0.058`$. This compares well with a single-giant KIII1-5 classification ($`R=1128R_{}`$ and $`\rho =134\times 10^4`$ g cm-3). Thus, the luminosity class - giant - indicates that the star’s envelope can be contained within the Roche lobe, and thus evolutionary expansion of the donor can sustain Roche lobe mass transfer rate. Wind-fed accretion alone cannot power the luminosities we observe from this system.

The long orbital period suggests that the evolution of GRS 1915+105 is further advanced than V404 Cyg, a similar system with a stripped-subgiant K0IV donor at a period of 6.5 days around the black hole ($`M_\mathrm{d}=0.7_{0.2}^{0.3}M_{}`$; Shahbaz et al. 1994). The more evolved giant star in GRS 1915+105 is also indicated by the smaller mass of the donor star compared to V404 Cyg due to increased mass loss sustained by nuclear evolution (King et al. 1996). Indeed, the mass of the donor star, a “stripped-giant”, is lower by a factor of three compared to a single giant star. King (1993) has shown that the mass of the donor star can have a maximum mass of $`1.3M_{}`$ while its actual mass depends on the mass of the helium core. Adapting the relation $`R_\mathrm{d}=12.55R_{}(M_\mathrm{c}/0.25)^{5.1}`$ from King (1993), we derive $`M_\mathrm{c}=0.27M_{}`$ for $`R_\mathrm{d}=19R_{}`$. The helium core has a mass $`M_\mathrm{c}=0.27M_{}`$, one third of the remaining mass of the star at $`M_\mathrm{d}=0.81M_{}`$. More than 1 $`M_{}`$ mass of the donor star has already been transferred to the black hole. The orbital period of GRS 1915+105 is increasing and will reach at least 45 days (section 3 in King 1993), a point when the donor star will reach its helium-core mass and the mass-transfer rate will cease. Actually, the stellar size for the helium-core mass of $`M_\mathrm{c}=0.27M_{}`$ is $`R_\mathrm{d}=18.6R_{}`$ (relation (1) in King 1993), suggesting that the stellar size is close to its possible minimum, thus not far from the time the mass transfer rate will cease.

GRS 1915+105 is in perpetual outburst since its discovery in 1992 whereas V404 Cyg has been in quiescence after its X-ray outburst. This is consistent with a higher mass transfer rate expected towards the end of the nuclear evolution of the donor star. Indeed, a higher mass-transfer rate is expected in GRS 1915+105 by 4–8 times compared to V404 Cyg, by using the relation, $`\dot{M}_25.4\times 10^9(M_\mathrm{c}/0.25)^{7.11}M_\mathrm{d}M_{}`$ yr-1 (King 1993), compared to an Eddington accretion-rate of $`\dot{M}_{\mathrm{Edd}}1.5\times 10^8`$ erg s-1 (for a 14 $`M_{}`$ black-hole mass). An increased orbital period, for example from the 6.5 days of V404 Cyg to the 33.5 days of GRS 1915+105, due to evolution would result in longer outbursts and a decreased duty cycle as well as higher super-Eddington accretion rates (King et al. 2002). The observed X-ray luminosity during outburst is dependent on inclination and anisotropic radiation patterns, therefore it is not easy to show that the accretion rate is higher in GRS 1915+105. The observed X-ray luminosity is super-Eddington at the outburst maximum ($`L_{\mathrm{Edd}}1.8\times 10^{39}`$ erg s-1 for a 14 $`M_{}`$ black hole mass) in both V404 Cyg at $`L_\mathrm{x}(\mathrm{max})3.75\times 10^{39}(D/5\mathrm{kpc})^2`$ erg s-1 (Tanaka 1989) and in GRS 1915+105 at $`L_\mathrm{x}(\mathrm{max})5\times 10^{39}(D/12\mathrm{kpc})^2`$ erg s-1 (Rau et al. 2003). Interestingly, King (2002) recently proposed that GRS 1915+105, a low-mass X-ray binary showing super-Eddington accretion, may well be the best representative in our galaxy of ultraluminous X-ray sources in elliptical galaxies.

It is clear that, as far as the physical mechanisms involved in the radiation emitted, this rich system will be explored more efficiently once the binary parameters are known more accurately. A better sampling of the orbital period will decrease the uncertainty in the stars’ masses which mostly arises from the radial velocity uncertainty, as well as the independent determination of the binary inclination from the IR ellipsoidal modulations for comparison to the inclination derived from the radio jets will constrain further evolutionary considerations of this peculiar system.

###### Acknowledgements.

We thank Bart Willems for discussion on tidal interaction in close binaries and an anonynous referee for helpful comments. The use of the MOLLY software package developed by T. Marsh is acknowledged.


text:
# e-print published at http://www.arxiv.org/abs/cond-mat/0312317 December 2003 Estimating Driving Forces of Nonstationary Time Series with Slow Feature Analysis

## 1 Introduction

Nonlinear time series analysis generally assumes stationarity (see Casdagli, , 1997; Schreiber, , 1999, for an overview). However, many time series are actually nonstationary for various reasons, such as temperature drift in the experimental setup, decreasing reservoirs in (bio)chemical reactors or ecological systems, global warming in climate data, varying heart rate in cardiology, or vibrato/tremor in vocal fold vibrations. Such nonstationarities can be modeled by underlying parameters, referred to as driving forces, that change the dynamics of the system smoothly on a slow time scale or abruptly but rarely, e.g. if the dynamics switches between different discrete states.

If a test reveals that a time series is nonstationary Witt et al.,  (1998), one can still apply methods of stationary time series analysis if one determines sections where the driving force has similar values and analyzes these sections as one stationary time series. This can be done by first slicing the time series into windows of equal size and then grouping the windows based on similarity measures of the dynamics Manuca & Savit,  (1996). This division of the time series can be avoided by the technique of overembedding. If the embedding dimension is sufficiently high, then similar embedding vectors automatically belong to similar values of the driving forces and all available data can be used for analysis, such as forecasting or nonlinear noise reduction Hegger et al.,  (2000).

However, in some cases, e.g. in the analysis of EEG data, one is particularly interested in revealing the driving forces themselves, which is in principle only possible up to an invertible transformation. One standard method for visualizing driving forces is the recurrence plot Eckmann et al.,  (1987), but it is often difficult to interpret. Methods have been developed to estimate driving forces based on the finding that the recurrence plot of a time series is similar to the recurrence plot of its underlying driving forces Casdagli,  (1997); Schreiber,  (1999). Another technique for the reconstruction of driving forces has been presented in Verdes et al.,  (2001). Here I present an alternative approach based on slow feature analysis, a new technique developed in the field of theoretical neurobiology.

## 2 Slow Feature Analysis

Slow feature analysis (SFA) has been originally developed in context of an abstract model of unsupervised learning of invariances in the visual system of vertebrates Wiskott,  (1998) and is described in detail in Wiskott & Sejnowski,  (2002).

The general objective of SFA is to extract slowly varying features from a quickly varying signal. For a scalar output signal it can be formalized as follows. Let $`𝒙=𝒙(t)`$ be an $`N`$-dimensional input signal where $`t`$ indicates time and $`𝒙=[x_1,\mathrm{},x_N]^T`$ is a vector. Find the input-output function $`g(𝒙)`$ that generates a scalar output signal

$$y(t):=g(𝒙(t))$$
(1)
with most slowly temporal variation possible as measured by the variance of the time derivative:

$$\text{minimize}\dot{y}^2$$
(2)
with $``$ indicating the temporal mean. For convenience and to avoid the trivial constant solution the output signal has to meet the following constraints:

$`y`$ $`=`$ $`0\text{(zero mean)},`$ (3)
$`y^2`$ $`=`$ $`1\text{(unit variance)}.`$ (4)
This is an optimization problem of variational calculus and as such difficult to solve. However, if we constrain the input-output function to be a linear combination of some fixed and possibly nonlinear basis functions, the problem becomes tractable and can be solved in the following way. Let $`𝒉^{}=𝒉^{}(𝒙)`$ be a vector of some fixed basis functions. To be concrete assume $`𝒉^{}`$ contains all monomials of degree 1 and 2. Applying $`𝒉^{}`$ to the input signal $`𝒙(t)`$ yields the nonlinearly expanded signal $`𝒛^{}(t)`$:

$`𝒉^{}(𝒙)`$ $`:=`$ $`[x_1,x_2,\mathrm{},x_N,x_1^2,x_1x_2,\mathrm{},x_N^2]^T,`$ (5)
$`𝒛^{}(t)`$ $`:=`$ $`𝒉^{}(𝒙(t)).`$ (6)
Assume $`g(𝒙)`$ is a linear combination of the basis functions plus a constant $`c`$, generating the output signal $`y(t)=g(𝒙(t))`$:

$`g(𝒙)`$ $`:=`$ $`c+𝒂_{}^{}{}_{}{}^{T}𝒉^{}(𝒙),`$ (7)
$`y(t)`$ $`=`$ $`c+𝒂_{}^{}{}_{}{}^{T}𝒛^{}(t).`$ (8)
Since $`𝒉^{}`$ contains all monomials of degree 1 and 2, $`g(𝒙)`$ can be any polynomial of degree 2, if the constant $`c`$ and the coefficient vector $`𝒂^{}`$ are chosen correspondingly.

For reasons that become clear below, it is convenient to sphere (or whiten) the expanded signal and transform the basis functions accordingly:

$`𝒉(𝒙)`$ $`:=`$ $`𝑺(𝒉^{}(𝒙)𝒛^{}),`$ (9)
$`𝒛(t)`$ $`:=`$ $`𝑺(𝒛^{}(t)𝒛^{}),`$ (10)
with the sphering matrix $`𝑺`$ chosen such that the signal components have a unit covariance matrix, i.e. $`𝒛𝒛^T=𝑰`$, which can be easily done with the help of principal component analysis (PCA). The signal components have also zero mean, $`𝒛=\mathrm{𝟎}`$, since the mean values have been subtracted.

The input-output function and the output signal can now be written in these transformed basis functions and sphered expanded signal, respectively:

$`g(𝒙)`$ $`=`$ $`𝒂^T𝒉(𝒙),`$ (11)
$`y(t)`$ $`=`$ $`𝒂^T𝒛(t).`$ (12)
This is a form in which the constraints (3, 4) are particularly easy to meet, since we find for any coefficient vector $`𝒂`$ with norm 1,

$`y`$ $`=`$ $`𝒂^T\underset{=\mathrm{𝟎}}{\underset{}{𝒛}}=0,`$ (13)
$`y^2`$ $`=`$ $`𝒂^T\underset{=𝑰}{\underset{}{𝒛𝒛^T}}𝒂=𝒂^T𝒂=1,`$ (14)
which also motivates why constant $`c`$ has been dropped in (11, 12).

Thus the optimization problem reduces to finding the normalized coefficient vector $`𝒂`$ that minimizes

$$\dot{y}^2=𝒂^T\dot{𝒛}\dot{𝒛}^T𝒂,$$
(15)
which is obviously the normalized eigenvector of the time-derivative covariance matrix $`\dot{𝒛}\dot{𝒛}^T`$ with the smallest eigenvalue, which again can be easily found by PCA. The output signal $`y(t)`$ is then given by (12). Notice that $`y(t)`$ is usually uniquely determined up to the sign.

In summary SFA consists of basically four steps: (i) expand the input signal with some set of fixed possibly nonlinear functions; (ii) sphere the expanded signal to obtain components with zero mean and unit covariance matrix; (iii) compute the time derivative of the sphered expanded signal and determine the normalized eigenvector of its covariance matrix with the smallest eigenvalue; (iv) project the sphered expanded signal onto this eigenvector to obtain the output signal. The input-output function is given by (11) but is not needed here, since we are only interested in the extracted output signal.

In practice one has to work with time series $`𝒙_i`$ instead of continuous signals $`𝒙(t)`$, but the transfer of the algorithm described above to time series is straight forward. The time derivative is simply computed as the difference between successive data points without wrap-around assuming a constant sampling spacing $`\mathrm{\Delta }t`$.

The description of SFA given above assumes that all monomials of degree 1 and 2 are used for signal expansion. However, any other set of basis functions could be used as well, e.g. monomials of higher degree or radial basis functions. It is also straight forward to extend SFA to more than one output signal component, in which case additional components should be uncorrelated to earlier ones. Under this constraint additional components can be computed by using the normalized eigenvectors of (15) with the next larger eigenvalues Wiskott & Sejnowski,  (2002). SFA requires that all components vary to some extent. Thus if there are constant components in the input signal or the expanded signal, these are simply discarded.

## 3 Examples

In the following I will present two examples with time series $`w_i`$ derived from a tent map and a logistic map to illustrate the properties of SFA. I have also done simulations with the Lorenz System, but results were very unreliable and sensitive to parameter variations and noise. The underlying driving force will always be denoted by $`\gamma `$ and may vary between $`1`$ and $`1`$ either smoothly or rarely, but with a comparable time scale of variation (as defined by the variance of its time derivative (2)).

Taking the time series $`w_i`$ directly as an input signal would not give SFA enough information to estimate the driving force, because SFA considers only data (and its derivative) from one time point at a time. Thus it is necessary to generate an embedding-vector time series as an input. Here embedding vectors at time point $`i`$ with delay $`\tau `$ and dimension $`m`$ are defined by

$`𝒙_i`$ $`:=`$ $`[w_{i\tau (m1)/2},w_{i\tau ((m1)/21)},\mathrm{},w_{i+\tau (m1)/2}]^T,`$ (16)
for scalar $`w_i`$ and odd $`m`$. The definition can be easily extended to even $`m`$, which requires an extra shift of the indices by $`\tau /2`$ or its next lower integer to center the used data points at $`i`$. Centering the embedding vectors results in an optimal temporal alignment between estimated and true driving force.

The following simulations are based on 6000 data points each and were done with Matlab (Release 13).

### 3.1 Tent Map

As a first example consider a time series generated by an iterated tent map Casdagli,  (1997)

$`\widehat{\gamma }_i`$ $`:=`$ $`(\gamma _i+1)/2`$ (17)
$`w_{i+1}`$ $`=`$ $`\{\begin{array}{cccc}+2w_i+\widehat{\gamma }_i\text{if}\hfill & \hfill 0& w_i& \frac{1}{2}(1\widehat{\gamma }_i)\hfill \\ 2w_i\widehat{\gamma }_i+2\text{if}\hfill & \hfill \frac{1}{2}(1\widehat{\gamma }_i)<& w_i& \frac{1}{2}(2\widehat{\gamma }_i)\hfill \\ +2w_i+\widehat{\gamma }_i2\text{if}\hfill & \hfill \frac{1}{2}(2\widehat{\gamma }_i)<& w_i& 1\hfill \end{array},`$ (21)
which maps the interval $`[0,1]`$ onto itself with a functional form “$``$” for $`\gamma =1`$ (or $`\widehat{\gamma }=0`$). The parameter $`\gamma `$ shifts this wedge cyclically to the left until it becomes a “$``$” for $`\gamma =\widehat{\gamma }=+1`$. Figure 1 shows the true driving force, the time series, and the estimated driving force with $`m=10`$, $`\tau =1`$, and third order polynomials for SFA. The correlation between true and estimated driving force is $`r=0.96`$. Note that the scale and offset of the estimated driving force are arbitrarily fixed by the constraints and that the sign is random. The axes were therefore chosen such that the curves in the bottom graphs of this and the following figures are optimally aligned with each other.

### 3.2 Logistic Map

As a second example consider a time series derived from a logistic map

$$w_{i+1}=(3.6+0.4\gamma _i)w_i(1w_i),$$
(22)
which maps the interval $`[0,1]`$ onto the interval $`[0,0.9+0.1\gamma _i]`$ and has the shape of an upside-down parabola crossing the abscissa at 0 and 1. Parameter $`\gamma `$ governs the height of the parabola. Figure 2 shows the true driving force, the time series, and the estimated driving force with $`m=3`$, $`\tau =1`$, and second order polynomials for SFA. The correlation between true and estimated driving force is $`r=0.997`$.

## 4 Technical Remarks

#### Choice of Parameters

SFA is basically parameter-free except for the general choice of the class of nonlinear basis functions. The only other parameters to choose in this method are the dimension $`m`$ and the time delay $`\tau `$ of the embedding vectors. The value of $`\dot{y}^2`$ is a fairly reliable indicator for a good choice of values for $`m`$ and $`\tau `$. The smaller $`\dot{y}^2`$ the better, since there is no trivial way of achieving small $`\dot{y}^2`$-values (except if the number of basis functions comes close to the number of data points). I typically test a certain range of $`\tau `$-values and successively increase the $`m`$-value until performance (measured in terms of $`\dot{y}^2`$ or $`r`$) is satisfactory.

#### Rarely Varying Driving Forces

If SFA is designed to extract slowly varying features from a signal, how about driving forces that apparently violate this assumption heavily, e.g. if they vary abruptly and jump between different discrete values? The definition of slowness given in (2) does actually not depend on the smoothness of the output signal, it may equally well be a signal that switches abruptly between different discrete values as long as the jumps occur rarely enough to lead to the same variance of the time derivative. Whether the order of the values can be recovered depends on the dynamics of the driving force. SFA will tend to order the values such that jumps occur between nearby values of the estimated driving force. Thus, if the original values are -0.5, 0, and +0.5 and there are many direct jumps between -0.5 and +0.5, SFA might map the three values onto -1, +1, and 0, respectively, thereby changing the order of the second and third value. Figures 3 and 4 show results for the two systems if a rarely instead of a slowly and smoothly varying driving force is used.

#### High-Dimensional Input Data

Probably the most severe limitation of SFA is the fact that the number of monomials (or any other basis functions) grows quickly with the dimensionality of the embedding vectors (curse of dimensionality). One therefore might run out of computer memory before a high enough $`m`$-value is reached. However, higher-dimensional problems could be dealt with in a hierarchical fashion by breaking the embedding vectors into smaller parts which are first analyzed separately and the results of which are then combined for a final analysis. Applying this hierarchical scheme to 65-dimensional input vectors has been demonstrated in Wiskott & Sejnowski,  (2002).

#### Accuracy of the Estimated Driving Force

It might be surprising that in the examples above SFA is able to estimate the driving forces with such an accuracy up to a factor and a constant offset, even though the estimation is undetermined up to any invertible transformation, not only scaling and shift. One reason for that is that the driving forces used are relatively slow already and could probably not be improved much by an invertible nonlinear transformation. However, even if that were not the case one might hope that in practice for the lower-dimensional embedding vectors ($`m=10`$ instead of, e.g., $`m=50`$) SFA has to find a relatively simple input-output function $`g(𝒙)`$, which is more likely to preserve the exact shape of the driving-force curve than a more complicated one. Further experiments are needed to verify this.

#### Noise Sensitivity

How sensitive is the method to noise? One might suspect that SFA is quite sensitive to noise, since eigenvectors of smallest eigenvalues are used. However, focusing on small eigenvalues is only done after sphering and taking the time derivative, so that small but quickly varying noise components typically have large eigenvalues and are discarded by SFA. Graceful degradation with noise was also found in Wiskott & Sejnowski,  (2002). Thus, using the eigenvectors with smallest eigenvalues in itself does not induce noise sensitivity. I have done simulations with the examples presented here by adding Gaussian white noise to the signals before embedding. Adding 10%, 20%, and 50% noise to the tent map time series reduced the correlation between true and estimated driving force from $`r=0.96`$ to about 0.94, 0.90, and 0.71, respectively; adding 2%, 5%, and 10% noise to the logistic map time series reduced the correlation from $`r=0.997`$ to about 0.97, 0.87, and 0.70, respectively. Thus we see, that in these examples the method is in fact fairly robust with respect to noise.

#### Multiple Driving Forces

SFA can be easily extended to the extraction of multidimensional output signals Wiskott & Sejnowski,  (2002), which could be used to estimate multiple driving forces. However, this might require longer time series and higher-dimensional embedding vectors. If the different driving forces are not clearly separated by different time scales, SFA might only be able to estimate them up to a linear mixing transformation. In that case independent component analysis (ICA) might be able to separate them, if they are statistically independent and not more than one is Gaussian (see Hyvärinen, , 1999, for an overview of ICA methods).

## 5 Conclusion

In this paper I have demonstrated that slow feature analysis (SFA) can be applied to the problem of estimating a driving force of a nonstationary time series. The estimates are fairly accurate except for a scaling factor and a constant offset, which cannot be extracted in general. The method works for slowly as well as rarely varying driving forces, although in the latter case the order of the different discrete values might not be recoverable. The slowness principle also provides robustness to noise, which is typically quickly varying and therefore suppressed by SFA as far as possible.

There are still a number of open questions. The next steps of investigation will have to include a comparison with standard methods, such as the technique of recurrence plots, and an exploration of the conditions under which SFA can be applied with similar success as demonstrated here. It is also necessary to apply SFA to real world data. This has been successfully done in the context of learning receptive field properties of the visual cortex based on natural image sequences Berkes & Wiskott,  (2003), but that was not for extracting driving forces. In any case, since SFA works very differently from other techniques that have been applied to the estimation of driving forces, one can hope that it at least complements these other techniques in some cases.

## 6 Acknowledgment

I am grateful to Hanspeter Herzel and Isao Tokuda for useful hints and Hanspeter Herzel also for critically reading the manuscript. This work has been supported by the VolkswagenStiftung.


text:
# Noisy Control, the Adiabatic Geometric Phase, and Destruction of the Efficiency of Geometric Quantum Computation

## Abstract

We examine the adiabatic dynamics of a quantum system coupled to a noisy classical control field. A stochastic phase shift is shown to arise in the off-diagonal elements of the system’s density matrix which can cause decoherence. We derive the condition for onset of decoherence, and identify the noise properties that drive decoherence. We show how this decoherence mechanism causes: (1) a dephasing of the observable consequences of the adiabatic geometric phase; and (2) the loss of computational efficiency of the Shor algorithm when run on a sufficiently noisy geometric quantum computer.

As the adiabatic geometric phase approaches its twentieth birthday bry ; gph  it continues to impact research in physics and chemistry. The original scenario considered a quantum system $`𝒬`$ with a discrete, nondegenerate energy spectrum whose dynamics was driven by a classical set of control parameters $`𝐁(t)`$. It was shown that when $`𝐁(t)`$ was cycled adiabatically through a loop in parameter space, and $`𝒬`$ was initially prepared in an eigenstate $`|E(0)`$ of the initial Hamiltonian $`H(0)`$, its quantum state would return to the initial state $`|E(0)`$ at the end of the cycle, to within a phase factor $`\mathrm{exp}[i\varphi ]`$. The phase shift $`\varphi `$ contained a geometric contribution now known as the adiabatic geometric phase (AGP) which had been discarded in prior treatments of quantum adiabatic dynamics. Since then, this original scenario has been generalized in a number of ways: (1) the requirement of nondegeneracy was removed wze ; (2) the adiabatic restriction was removed AhA ; (3) the control field $`𝐁(t)`$ was allowed to have its own dynamics dyn ; and (4) $`𝒬`$ was allowed to interact with a quantum environment as well as with the control field $`𝐁(t)`$ qe1 ; qe2 ; qe3 . In this Letter we return to the original scenario, though we allow the classical control field to contain a noise component $`𝐁_n(t)`$ along with the (deterministic) adiabatically varying component $`𝐁_a(t)`$: $`𝐁(t)=𝐁_a(t)+𝐁_n(t)`$. We will show that $`𝐁_n(t)`$ causes a stochastic phase shift to appear in the off-diagonal elements of $`𝒬`$’s density matrix which, for sufficiently strong noise, causes decoherence. We derive the condition that determines the onset of decoherence (see eq. (6)), and identify the noise properties that drive the decoherence. We show that this decoherence mechanism causes all physical consequences of the AGP to become unobservable.

The AGP has also had an impact on the new research area of quantum computation. In 2000, a proposal was put forth for a geometric quantum computer (GQC) gqc  which uses the AGP to encode conditional phase shifts into the quantum state of a GQC. It was argued that the AGP gave the GQC a modest degree of fault tolerance to uncertainty in the control parameters, though the need to examine the robustness of a GQC to decoherence was also noted. In the second part of this Letter we carry out such a decoherence analysis. We show that the above decoherence mechanism arising from noisy control severely impacts the performance of a quantum algorithm run on a noisy GQC. Specifically, we show how noisy control causes the Shor algorithm to loose its computational efficiency so that the time to factor an integer ($`N`$) grows exponentially with the size of the integer ($`\mathrm{log}N`$). We establish the condition for the onset of decoherence, and show that the success probability for the algorithm is smoothly degraded as the noise becomes stronger, reaching the point at onset where the computational efficiency of the Shor algorithm has been completely destroyed. We also show that at onset of decoherence, the entanglement fidelity of a 2-qubit entangled state is degraded to that of a mixed state upon output from a sufficiently noisy controlled phase gate of the type used in a GQC. This noise-induced loss of entanglement confirms that the GQC is no longer behaving as a quantum computer. To the best of our knowledge, this Letter provides the first quantitative demonstration of how the computational efficiency of a landmark quantum algorithm such as Shor’s is destroyed when it is run on a sufficiently noisy GQC. A detailed presentation of all arguments will be reported elsewhere fg3 .

1. General Analysis: We consider the situation where the dynamics of an N-qubit system $`𝒬`$ is driven by an external control field $`𝐁(t)`$. The control field is composed of a deterministic component $`𝐁_a(t)`$ and a noise component $`𝐁_n(t)`$. We will be interested in cases where $`𝐁_a(t)`$ executes a cyclic evolution with period $`T`$. We take the noise component $`𝐁_n(t)`$ to be a stationary stochastic process with zero mean ($`\overline{𝐁_n(t)}=0`$) and noise correlation time $`\tau _c`$. To simplify the notation we assume the noise is isotropic. This assumption is not essential, and the following analysis can easily be repeated without it. For isotropic noise, the noise correlation function $`\overline{B_n^i(t)B_n^j(t^{})}=\delta _{ij}\sigma ^2f(\tau )`$. Here: (1) $`\delta _{ij}`$ is the Kronecker delta; (2) $`\sigma ^2`$ is the variance of the noise $`\overline{(B_n^i(t))^2}`$; (3) $`\tau =tt^{}`$; and (4) $`f(\tau )`$ is the normalized noise fluctuation profile with peak value $`f(0)=1`$ and temporal width of order $`\tau _c`$. Because of the stochastic nature of $`𝐁_n(t)`$, each application $`k`$ of the control field contains a different realization of the noise $`𝐁_n(t;k)`$. To account for the dynamical effects of noise, it is necessary to introduce an ensemble of identical quantum systems $`𝒬_k`$, all prepared in the same initial state $`|\psi (0)`$, but whose evolution is subject to a different noise realization $`𝐁_n(t;k)`$.

Quite generally, the system Hamiltonian $`H(t)`$ can be written as the sum of two terms: $`H_a(t)`$ and $`H_s(t)`$. Here $`H_a(t)`$ is the Hamiltonian in the absence of noise and depends solely on $`𝐁_a(t)`$. The stochastic interaction term $`H_s(t)`$ depends explicitly on the noise component $`𝐁_n(t)`$ and its action varies from one application of the control field to another. We introduce the instantaneous energy eigenvalues $`E_n(t)`$ and eigenstates $`|E_n(t)`$ of the noiseless Hamiltonian $`H_a(t)`$: $`H_a(t)|E_n(t)=E_n(t)|E_n(t)`$. Note that the periodicity of $`𝐁_a(t)`$ implies that $`H_a(t+mT)=H_a(t)`$. As in the original AGP scenario, the instantaneous energies $`E_n(t)`$ are assumed to be discrete and nondegenerate for all $`t`$ of interest so that $`|E_n(t+mT)=|E_n(t)`$. Introducing the minimum energy level separation $`\mathrm{\Delta }=\mathrm{min}|E_n(t)E_m(t)|`$, the noiseless dynamics will be adiabatic if $`\mathrm{}/T\mathrm{\Delta }`$. Furthermore, as shown in Ref. qe1 , the dynamics associated with the noise will also be adiabatic if $`\mathrm{}/\tau _c\mathrm{\Delta }`$. Both of these conditions are assumed in the remainder of this Letter.

As noted above, each ensemble element $`𝒬_k`$ sees a different noise realization $`𝐁_n(t;k)`$ added to $`𝐁_a(t)`$. The dynamics of $`𝒬_k`$ is driven by $`H(t;k)=H_a(t)+H_s(t;k)`$ which generates the final state $`|\psi (t;k)=U(t_f,t_0;k)|\psi (0)`$, where $`U(t_f,t_0;k)=\mathrm{exp}\left[(i/\mathrm{})_{t_0}^{t_f}𝑑tH(t;k)\right]`$. Thus the final state produced by the noisy control field varies over the ensemble. Consequently, calculation of the expectation value of an observable requires an ensemble average over the noise, along with the usual quantum mechanical averaging.

We divide the time interval ($`t_0`$, $`t_f`$) into $`J`$ subintervals of duration $`ϵ=(t_ft_0)/J`$ by introducing intermediate times $`t_j=t_0+jϵ`$ ($`j=0,\mathrm{},J`$). Eventually we let $`J\mathrm{}`$. The propagator $`U(t_f,t_0)`$ factors into a product of propagators: $`U(t_f,t_0)=U(J)\mathrm{}U(0)`$, where $`U(j)=U(t_j,t_{j1})`$. One can show fg3 ; fg4  that in the adiabatic limit the matrix elements $`U_{lm}(j)=E_l(t_j)|U(j)|E_m(t_{j1})`$ are:

$`U_{lm}(j)=`$ (1)
$`\delta _{lm}\left[1{\displaystyle \frac{iϵ}{\mathrm{}}}\left\{E_l(t_j)\mathrm{}\dot{\gamma }_l(t_j)+H_s^{ll}(t_j)\right\}+𝒪(ϵ^2)\right].`$
Here $`\dot{\gamma }_l(t_j)=iE_l(t_j)|\dot{E}_l(t_j)`$ is the time derivative of the AGP for energy level $`l`$; $`\mathrm{\Gamma }_{lm}(t_j)=E_l(t_j)|\dot{E}_m(t_j)`$ is the non-adiabatic coupling for levels $`l`$ and $`m`$; and $`H_s^{ll}(t_j)=E_l(t_j)|H_s(t_j)|E_l(t_j)`$. Choosing $`|\psi (t_0)=|E_k(t_0)`$; reconstructing $`U(t_f,t_0)`$ from the $`U(j)`$; and letting $`J\mathrm{}`$ gives:

$$|\psi (t_f)=\mathrm{exp}\left[i\left\{\mathrm{\Gamma }_a(k)+\mathrm{\Gamma }_s(k)\right\}\right]|E_k(t_0).$$
(2)
Here $`\mathrm{\Gamma }_a(k)=(1/\mathrm{})_{t_0}^{t_f}𝑑t\left\{E_k(t)\mathrm{}\dot{\gamma }_k(t)\right\}`$ and $`\mathrm{\Gamma }_s(k)=(1/\mathrm{})_{t_0}^{t_f}𝑑tH_s^{kk}(t)`$. Eq. (2) shows that noise introduces a stochastic phase shift $`\mathrm{\Gamma }_s(k)`$ into the final state $`|\psi (t_f)`$ which is sensitive to the initial energy eigenstate $`|E_k(t_0)`$ through $`H_s^{kk}(t)`$.

To bring out the physical consequences of the stochastic phase $`\mathrm{\Gamma }_s(k)`$ we consider an initial superposition state $`|\psi (t_0)=_kc_k|E_k(t_0)`$. From the linearity of quantum mechanics: $`|\psi (t_f)=_kc_kU(t_f,t_0)|E_k(t_0)`$. Information about the phase coherence of the final superposition is carried in the off-diagonal elements of the final density operator $`\rho (t_f)=|\psi (t_f)\psi (t_f)|`$. One can show fg3  that the noise averaged density matrix elements are:

$$\overline{\rho _{kj}(t_f)}=c_kc_j^{}\mathrm{exp}\left[i\mathrm{\Gamma }_a(k,j)\right]\overline{D(k,j)},$$
(3)
where $`\mathrm{\Gamma }_a(k,j)=\mathrm{\Gamma }_a(k)\mathrm{\Gamma }_a(j)`$, and the decoherence factor $`\overline{D(k,j)}`$ is:

$$\overline{D(k,j)}=\overline{\mathrm{exp}\left[i\mathrm{\Gamma }_s(k,j)\right]},$$
(4)
with $`\mathrm{\Gamma }_s(k,j)=\mathrm{\Gamma }_s(k)\mathrm{\Gamma }_s(j)`$. We see that the decoherence factor $`\overline{D(k,j)}`$ is a direct consequence of the noise-induced phase shift $`\mathrm{\Gamma }_s`$. Often, $`H_s(t)`$ is linear in the noise component $`𝐁_n(t)`$ so that $`H_s(t)=(\gamma \mathrm{}/2)𝐁_n(t)\widehat{𝐎}`$. Here $`\gamma `$ is the coupling constant and $`\widehat{𝐎}`$ is a vector operator. In the adiabatic limit, $`\tau _c(t_ft_0)𝒯`$. Thus, if we partition the integration interval ($`t_0`$, $`t_f`$) that appears in the definition of $`\mathrm{\Gamma }_s(k,j)`$ into $`M=𝒯/\tau _c`$ subintervals of duration $`\tau _c`$, we render $`\mathrm{\Gamma }_s(k,j)`$ into a sum of uncorrelated random variables $`\mathrm{\Gamma }_s^m(k,j)`$: $`\mathrm{\Gamma }_s(k,j)=_{m=1}^M\mathrm{\Gamma }_s^m(k,j)`$. Because the noise is stationary, the $`\mathrm{\Gamma }_s^m(k,j)`$ have identical probability distributions. If the $`\mathrm{\Gamma }_s^m(k,j)`$ are not only uncorrelated, but also statistically independent, it follows from the Central Limit Theorem that $`\mathrm{\Gamma }_s(k,j)`$ will have a Gaussian probability distribution. In this case, one can show fg3  that $`\overline{\mathrm{\Gamma }_s(k,j)}=0`$ and

$$\overline{\mathrm{\Gamma }_s^2(k,j)}=\frac{\eta \gamma ^2\sigma ^2}{4}I_{kj},$$
(5)
where: (1) $`\eta `$ cycles of the control field have been applied ($`\eta T=𝒯`$); (2) $`I_{kj}=_0^T𝑑t𝑑t^{}\widehat{𝐎}_{kj}(t)\widehat{𝐎}_{kj}(t^{})f(tt^{})`$, with $`\widehat{𝐎}_{kj}(t)=\widehat{𝐎}_{kk}(t)\widehat{𝐎}_{jj}(t)`$, and $`\widehat{𝐎}_{kk}(t)=E_k(t)|\widehat{𝐎}|E_k(t)`$; and (3) $`f(tt^{})`$ is the normalized noise fluctuation profile introduced earlier. Having the mean and variance of the Gaussian probability distribution for $`\mathrm{\Gamma }_s(k,j)`$, we can evaluate the noise average in eq. (4). The result is $`\overline{D(k,j)}=\mathrm{exp}[\overline{\mathrm{\Gamma }_s^2(k,j)}/2]`$. One expects that when the phase uncertainty $`\sqrt{\overline{\mathrm{\Gamma }_s^2(k,j)}}2\pi `$, the superposition in $`|\psi (t_f)`$ will have been dephased/decohered by the noise. Inserting $`\overline{\mathrm{\Gamma }_s^2(k,j)}(2\pi )^2`$ into our expression for $`\overline{D(k,j)}`$ gives $`\overline{D(k,j)}3\times 10^9`$ so that the off-diagonal elements of $`\overline{\rho (t_f)}`$ (see eq. (3)) are effectively zero and the noise has in fact caused an effective collapse of the wavefunction $`|\psi (t_f)`$. One can show fg3  that the noise variance $`\sigma ^2`$, the average noise power absorbed by $`𝒬`$ per unit volume $`\overline{P}/V`$, and the effective bandwidth $`\mathrm{\Delta }\omega `$ of the absorbed noise power are related: $`\overline{P}/V=\sigma ^2\mathrm{\Delta }\omega `$. Using this relation, together with eq. (5), one can re-express the condition for onset of decoherence $`\sqrt{\overline{\mathrm{\Gamma }_s^2(k,j)}}2\pi `$ as

$$\frac{\eta }{16\pi ^2}\left(\frac{\gamma ^2}{\mathrm{\Delta }\omega }\right)\left(\frac{\overline{P}}{V}\right)I_{jk}1.$$
(6)
We consider a number of applications of this general analysis in the remainder of this Letter.

2. Dephasing the AGP: In our first application we reconsider the original AGP scenario bry , allowing for noisy control. Here Q is a single qubit, and $`𝐁_a(t)`$ precesses about the z-axis with period $`T`$ at an angle $`\theta _0`$. Let $`C`$ denote the contour traced out by the tip of $`\widehat{𝐁}_a(t)`$ in a time $`T`$. To observe the AGP, the initial state must be a superposition: $`|\psi (0)=c_+|E_+(0)+c_{}|E_{}(0)`$. Eq. (3) gives the matrix elements of the final density matrix. The diagonal elements are real-valued, and thus contain no information about the AGP. The off-diagonal elements, however, do depend on the AGP through $`\mathrm{\Gamma }_a(+,)=\mathrm{\Gamma }_a(+)\mathrm{\Gamma }_a()`$, where $`\mathrm{\Gamma }_a(\pm )`$ are defined below eq. (2). In an NMR experiment, observation of the AGP is carried out by measuring the transverse magnetization whose expectation value depends on $`\overline{\rho _+(t_f)}`$ sut ; fg4 . As the analysis of Section 1 showed, sufficiently noisy control causes $`\overline{\rho _+(t_f)}`$ to effectively vanish so that all physical consequences of the AGP thus become unobservable. The stochastic phase shifts $`\mathrm{\Gamma }_s(\pm )`$ generated by noisy control dephase the final superposition state, reducing it to a mixture of the states $`|E_\pm (0)`$ in which all AGP effects are absent. In the NMR setting, the precession of $`𝐁_a(t)`$ is produced by varying the phase of the rf magnetic field $`𝐁_{rf}(t)=B_{rf}(t)\widehat{𝐱}`$. If we imagine the noise is due to rf power fluctuations, $`H_s(t)=(\gamma \mathrm{}/2)B_n(t)\sigma _x`$ which identifies $`\widehat{𝐎}=\widehat{𝐱}\sigma _x`$. One can show fg3  that for this type of noise $`I_+=4\tau _cT\mathrm{sin}^2\theta _0`$. If we denote the static magnetic field that splits the nuclear energy levels by $`𝐁_0=B_0\widehat{𝐳}`$, it is well-known sut  that $`\mathrm{sin}^2\theta _0=B_{rf}^2/(B_{rf}^2+(B_0B_{rf})^2)`$. Using this result for $`I_+`$ in eq. (6) gives:

$`{\displaystyle \frac{\eta \gamma ^2T}{4\pi ^2}}\left[{\displaystyle \frac{B_{rf}^2}{B_{rf}^2+(B_0B_{rf})^2}}\right]\left({\displaystyle \frac{\overline{P}}{V}}\right)\left({\displaystyle \frac{\tau _c}{\mathrm{\Delta }\omega }}\right)`$ $``$ $`1.`$ (7)
Eq. (7) identifies the noise properties that impact decoherence, and allows a quantitative assessment of how much noise can be tolerated before all AGP effects are dephased by noise. We believe that this is the first demonstration of how noisy control causes a dephasing of the AGP, and the first quantitative analysis of the onset condition for this dephasing mechanism.

3. Shor Algorithm on a Noisy GQC: We now consider how noisy control impacts the performance of a GQC. The universal set $`U`$ of quantum gates used to construct a GQC contains the 1-qubit Hadamard gate $`H`$ and the set of all possible 2-qubit controlled-phase gates $`B(\varphi )`$ ($`\varphi [0,2\pi )`$). The action of $`B(\varphi )`$ on the 2-qubit computational basis states (CBS) $`|xy`$ is: $`B(\varphi )|xy=\mathrm{exp}[ixy\varphi ]|xy`$, where x, y = 0, 1. The AGP is used to encode the conditional phase shift $`\varphi `$. We focus on the operation of $`B(\varphi )`$ in the presence of noise as it is the only gate in the set $`U`$ that can introduce entanglement into the dynamics of a GQC. The conditional phase shift $`\varphi `$ is implemented using a four part pulse sequence $`P=P_0P_1P_2P_3`$, with $`P_0`$ ($`P_3`$) applied first (last). Here $`P_0=C\pi _1`$; $`P_1=\overline{C}\pi _2`$; $`P_2=P_0`$; and $`P_3=P_1`$, where $`C`$ is the cyclic evolution of $`𝐁_a(t)`$ introduced in Section 2; $`\overline{C}`$ is the time-reverse of $`C`$; and $`\pi _1`$ ($`\pi _2`$) is a $`\pi `$-pulse applied to the first (second) qubit. By appropriate choice of $`C`$, any phase shift $`\varphi `$ can be produced. The 1-qubit CBS are: $`|0=|E_{}(t)`$ and $`|1=|E_+(t)`$. Let $`k=(i_1,i_2)`$ with $`i_1`$, $`i_2`$ = 0, 1; then $`|E_k(t_0)|E_{i_1}(t_0)|E_{i_2}(t_0)`$. The pulse sequence $`P`$ maps $`k\stackrel{P_0}{}k_1\stackrel{P_1}{}k_2\stackrel{P_2}{}k_3\stackrel{P_3}{}k_4`$, with $`k_j=(i_1_{l=1}^jlmod2,i_2_{l=0}^{j1}lmod2)`$, and $``$ is addition modulo 2. For $`|\psi (t_0)=|E_k(t_0)`$, eq. (2) gives the final state $`|\psi (t_f)`$ with $`\mathrm{\Gamma }_a(k)=(1/\mathrm{})_{l=0}^3\mathrm{\Gamma }_a(k_l)`$; $`\mathrm{\Gamma }_s(k)=(1/\mathrm{})_{l=0}^3\mathrm{\Gamma }_s(k_l)`$; and $`\mathrm{\Gamma }_a(k_l)=_{t_0(l)}^{t_f(l)}𝑑t[E_{k_l}(t)\mathrm{}\dot{\gamma }_{k_l}(t)]`$; $`\mathrm{\Gamma }_s(k_l)=_{t_0(l)}^{t_f(l)}𝑑tH_s^{k_lk_l}(t)`$, with $`t_0(l)=t_0+lT`$ and $`t_f(l)=t_0(l)+T`$. For the initial superposition state $`|\psi (t_0)=c_k|E_k(t_0)+c_j|E_j(t_0)`$, the final density matrix $`\overline{\rho _{kj}(t_f)}`$ is given by eq. (3). The decoherence factor $`\overline{D(k,j)}`$ is given by eq. (4) with $`\mathrm{\Gamma }_s(k,j)=_{l=0}^3\mathrm{\Gamma }_s^l(k,j)`$, and $`\mathrm{\Gamma }_s^l(k,j)=_{t_0(l)}^{t_f(l)}𝑑t[H_s^{k_lk_l}(t)H_s^{j_lj_l}(t)]`$. Because $`\tau _cT`$ in the adiabatic limit, the $`\mathrm{\Gamma }_s^l(k,j)`$ are uncorrelated so that $`\overline{\mathrm{\Gamma }_s^2(k,j)}=_{l=0}^3\overline{(\mathrm{\Gamma }_s^l(k,j))^2}`$. If the $`\mathrm{\Gamma }_s^l(k,j)`$ are also statistically independent, one can show fg3  that $`\overline{D(k,j)}=\mathrm{exp}[\overline{\mathrm{\Gamma }_s^2(k,j)}/2]`$ and $`\overline{\mathrm{\Gamma }_s^2(k,j)}`$ is given by eq. (5) with $`I_{kj}_{l=0}^3I_{kj}^l`$. Here $`I_{kj}^l`$ has the same form as $`I_{kj}`$, except that $`\widehat{𝐎}_{kj}(t)\widehat{𝐎}_{kj}^l(t)=\widehat{𝐎}_{k_lk_l}(t)\widehat{𝐎}_{j_lj_l}(t)`$. The condition for onset of decoherence is again $`\overline{\mathrm{\Gamma }_s^2(k,j)}4\pi ^2`$. As an example, imagine we input the Bell state $`(1/\sqrt{2})\left[|00+|11\right]`$ into $`B(\varphi )`$. Then $`c_k=c_j=1/\sqrt{2}`$; $`k=(00)`$; and $`j=(11)`$. If we assume again that the noise is due to rf power fluctuations, one can show fg3  that the noise-averaged entanglement fidelity $`\overline{F}=\overline{\psi (t_0)|\rho (t_f)|\psi (t_0)}=\frac{1}{2}+\frac{1}{2}\mathrm{cos}\mathrm{\Gamma }_a(k,j)\overline{D(k,j)}`$, with $`\mathrm{\Gamma }_a(k,j)=\mathrm{\Gamma }_a(k)\mathrm{\Gamma }_a(j)`$. Since $`\overline{D(k,j)}=\mathrm{exp}[\overline{\mathrm{\Gamma }_s^2(k,j)}/2]`$, eq. (5) shows that it goes smoothly to zero with increasing noise variance $`\sigma ^2`$. Thus $`\overline{F}1/2`$ which is the entanglement fidelity for an initial state which is a uniform mixture of the states $`|00`$ and $`|11`$, indicating that entanglement is destroyed once the control field driving $`B(\varphi )`$ becomes sufficiently noisy. One can show fg3  that for this Bell state $`_{l=0}^3I_{kj}^l=32\tau _cT\mathrm{sin}^2\theta _0`$ so that the condition for onset of decoherence (and loss of entanglement) is $`(2\eta /\pi ^2)(\gamma ^2/\mathrm{\Delta }\omega )(\overline{P}/V)(\tau _cT\mathrm{sin}^2\theta _0)1`$. Recent work t+n  also found that noise will severely impact the performance of $`B(\varphi )`$.

We now show how a GQC containing noisy controlled-phase gates impacts the computational efficiency of the Shor algorithm for factoring an integer $`N`$ shr . Number theoretic arguments reduce factoring to finding the period $`r`$ of the function $`F(a)=y^amodN`$, where $`y`$ is co-prime with $`N`$. The algorithm begins by preparing the GQC in the state $`|\psi _0=(1/\sqrt{A+1})_{j=0}^A|jr+l`$. Here $`r`$ is the period of $`F(a)`$; $`l`$ is the result of a measurement carried out during preparation of $`|\psi _0`$; and $`A`$ is the largest integer such that $`Ar+l<q`$, where $`q`$ is chosen such that $`N^2q2N^2`$. The algorithm implements a discrete Fourier transform modulo $`q`$ ($`DFT_q`$) on $`|\psi _0`$: $`|\psi _1=DFT_q|\psi _0=_{c=0}^{q1}\stackrel{~}{f}(c)|c`$, where

$$\stackrel{~}{f}(c)=\frac{\sqrt{r}}{q}\underset{j=0}{\overset{q/r1}{}}\mathrm{exp}\left[\frac{2\pi i}{q}(jr+l)c+i\mathrm{\Gamma }_s(c)\right].$$
(8)
The $`DFT_q`$ requires the application of $`L(L1)/2`$ controlled-phase gates, where $`L\mathrm{log}_2q`$. Noisy control of this gate causes the stochastic phase shift $`\mathrm{\Gamma }_s(c)`$ to appear in eq. (8). One can show fg3  that $`\overline{\mathrm{\Gamma }_s(c)}=0`$, and the variance $`\overline{\mathrm{\Gamma }_s^2(c)}`$ is given by eq. (5) with $`\eta =L(L1)/2`$; $`\sigma ^2=\overline{P}/V\mathrm{\Delta }\omega `$; and $`I_{kj}_{l=0}^3I_{kj}^l`$. The final step in the algorithm is a measurement which produces the result $`c`$ with noise-averaged probability $`P(c)=\overline{|\stackrel{~}{f}(c)|^2}`$. In the absence of noise, constructive interference in eq. (8) occurs when $`|rcc^{}q|r/2`$, which determines a unique $`c^{}`$, and consequently, $`P(c^{})=P(c)`$. The measurement result $`c`$ yields the period $`r`$ only if $`c^{}`$ and $`r`$ are co-prime. The success probability for the algorithm is then $`P_{suc}=_c^{}^{}P(c^{})`$, where in the primed sum only those $`c^{}`$ appear that are both less than, and co-prime with, $`r`$. One can show fg3  that:

$`\overline{|\stackrel{~}{f}(c)|^2}={\displaystyle \frac{1}{q}}`$ (9)
$`+{\displaystyle \frac{2r}{q^2}}{\displaystyle \underset{k=0}{\overset{q/r1}{}}}{\displaystyle \underset{j>k}{}}\mathrm{cos}\left[{\displaystyle \frac{2\pi (jk)}{q}}(rcmodq)\right]\overline{D(k,j)},`$
where $`\overline{D(k,j)}`$ is the decoherence factor due to noise. As above, $`\overline{D(k,j)}=\mathrm{exp}\left[\overline{\mathrm{\Gamma }_s^2(k,j)}/2\right]`$ which, by eq. (5), goes smoothly to zero as the noise variance $`\sigma ^2`$ increases. In the decoherence limit only the first term on the RHS of eq. (9) survives so that $`\overline{|\stackrel{~}{f}(c)|^2}1/q1/N^2`$. This result is independent of $`c`$ so that $`P_{suc}(1/N^2)\varphi (r)`$, where $`\varphi (r)`$ is Euler’s Phi function which gives the number of integers that are less than and co-prime with $`r`$. For large $`N`$, $`rN`$ and $`\varphi (r)r/\mathrm{log}rN/\mathrm{log}N`$. Thus, in the decoherence limit, $`P_{suc}1/(N\mathrm{log}N)=2^{\mathrm{log}N}/\mathrm{log}N`$. We see that the algorithm must be run on average $`2^{\mathrm{log}N}`$ times to obtain $`P_{suc}1`$. This is exponential in the problem size $`\mathrm{log}N`$, indicating that sufficiently strong noise destroys the computational efficiency of the Shor algorithm. This contrasts with the (computationally efficient) noiseless Shor algorithm whose runtime scales linearly with $`\mathrm{log}N`$. Onset of decoherence occurs when $`\overline{\mathrm{\Gamma }_s^2(c)}4\pi ^2`$, and for a GQC, $`_{l=0}^3I_{kj}^l\tau _cT\mathrm{sin}^2\theta _0`$ so that decoherence occurs when:

$$\left(\frac{\overline{P}}{V}\right)\left(\frac{\tau _c}{\mathrm{\Delta }\omega }\right)\frac{\pi ^2}{TL(L1)\gamma ^2\mathrm{sin}^2\theta _0}.$$
(10)
Adiabatic operation of a GQC requires $`T`$ to be large; and the desire to factor large $`N`$ means $`L=\mathrm{log}_2q\mathrm{log}_2N`$ will also be large. Thus the RHS of eq. (10) is expected to be small, so that managing noise will in fact be a significant issue for a GQC after all. Eq. (10) allows a quantitative estimate of how much noise can be tolerated by a GQC before decoherence due to noisy control undermines its operation. We believe that the above analysis is the first to quantitatively demonstrate how noisy control destroys the efficiency of Shor algorithm when run on a sufficiently noisy GQC.

###### Acknowledgements.

I would like to thank T. Howell III for continued support; the National Science Foundation for support through grant NSF-PHY-0112335; and the Army Research Office for support through grant DAAD19-02-1-0051.


text:
# Klio: A 5 micron camera for the detection of giant exoplanets

## 1 Why a 5 micron coronograph?

With the number of known extrasolar planets now exceeding 100 (http://www.obspm.fr/ encycl/catalog.html), we have reached the point where we can begin to study the statistical properties of these planets (i.e. Marcy et al. (2003)). However, while radial velocity surveys have provided a wealth of information about exoplanets, they are unable to determine the inclination of the systems or any spectral information. In addition, they are limited to relatively small star-planet separations over a practical baseline. Other indirect methods such as astrometry and microlensing are also valuable tools for the detection of exoplanets. However, only by direct detection can any spectral information about the exoplanets be determined, which is crucial to studying their physical properties (e.g. temperature, surface gravity, composition).

So far, the transit method is the only direct technique that has successfully detected an exoplanet Charbonneau et al. (2000) Konacki et al. (2003). While transit techniques clearly have the potential to provide a wealth of information about exoplanets, they are also limited to close in planets, where multiple transits can quickly be observed and exoplanets are more likely to transit their primary star.

In an effort to directly detect exoplanets at separations left largely unexplored by current measurements, we are building a coronographic imager to be used in the thermal infrared with ground-based adaptive optics (AO). The camera will be optimized to target the thermal IR (3-5 microns), where exoplanets are thought to peak in brightnessSudarsky, Burrows, & Hubeny (2003). Until now, it has been impossible to use the thermal IR since no appropriate space-based telescope exists and the sky background is so large from the ground. Additional emissivity introduced by conventional AO systems overwhelms any astronomical signal.

## 2 The adaptive secondary at the MMT

The adaptive secondary mirror at the 6.5m MMT provides a solution to the problem of increased emissivity from conventional adaptive optics systems. By making the secondary mirror of the telescope itself the deformable surface, approximately eight additional surfaces are eliminated compared to typical AO systems Lloyd-Hart (2000). This results in both an emissivity and throughput that are similar to a non-AO equipped telescope. For faint objects, a conventional AO system at the MMT would need to integrate 2-3 times longer than the adaptive secondary to achieve the same signal-to-noise in L and M bands Lloyd-Hart (2000).

The adaptive secondary mirror is now operational at the MMT telescope Wildi, Brusa, & Lloyd-Hart (2003) Brusa-Zappellini et al. (2003). Measurements taken with the BLINC/MIRAC mid-IR camera give an emissivity of 6.5% of the telescope plus adaptive optics system, which is typical of a telescope with no adaptive optics. In January 2003, BLINC/MIRAC also recorded the first high order AO images taken in the thermal IR. These were taken at 10.3 microns and show an achieve Strehl ratio (SR) of 0.96 with the AO system operational. Once vibration-related issues are resolved, the adaptive secondary should be able to achieve SR=0.9 at 5 microns. This number is based on actual observations at H band with the vibrations removed in post-processing.

## 3 The 5 micron coronograph

Klio will be a three channel system that consists of f/20, f/35, and pupil imaging modes. The primary design driver for Klio is use in the high background regime (i.e. L and M bands). The f/20 channel has been optimized for the L and M bands with Nyquist sampling at L band (0.048”/pixel) and a field of view of 15 x 12”. The f/35 channel is optimized for the H and K bands with Nyquist sampling at K band (0.027”/pixel) and a FOV of 8.7 x 7”.

Figure 1 shows the optical design of the f/20 channel. At the secondary focus of the telescope, a wheel will allow a selection of coronographic stops to be moved into the beam. Optics will then form an image of the secondary, where another wheel will allow a selection of Lyot stops to be introduced into the beam before the final image is formed

on the detector. The detector is an Indigo 320x256 pixel InSb chip with 30 micron pixels. We plan to have an engineering run without AO in the summer of 2004 and an engineering run with AO in the winter of 2004.

## 4 Exoplanet survey

So far, the only spectral information available for an exoplanet comes from observations of HD 209458b Charbonneau et al. (2001). As a result, our understanding of exoplanets, especially those at Jupiter-like separations, would benefit greatly from even a single direct detection. With ages derived from their primary stars, IR fluxes would provide anchor points for models that are essentially untested. In addition, a carefully constructed survey has the ability to provide statistical information about the exoplanet population, even for a null result.

Radial velocity surveys provide a baseline from which we can make predictions about the distribution of exoplanets at further separations from their primaries. Currently, results indicate a companion mass distribution of $`dn/dMM^{0.7}`$ and a separation distribution of $`dn/dlogaa^{>0}`$ Marcy et al. (2003). A fit to the data (for 0.2 AU ¡ a ¡ 2 AU) gives $`dn/daa^{0.7}`$ as a starting point. Our survey will enable a quantitative test of the hypothesis that the mass and separation distribution of outer planets are consistent with those observed within 5 AU.

As a first order prediction, we can simply extend the observed distributions to larger separations. We performed a series of Monte Carlo simulations where the radial velocity mass distribution was taken to be valid from 1-15 $`M_{Jupiter}`$ and the separation distribution was extended from 0.01-50 AU. An outer boundary of 50 AU was chosen to be consistent with the outer edge of our Kuiper Belt Trujillo & Brown (2001) as well as the radius of the smallest silhouette disk in Orion Bally, O’Dell, & McCaughrean (2000). Estimates of our sensitivities versus separation for L and M band are given in Heinze, Hinz, & McCarthy (2003) Heinze, Hinz, and McCarthy (2003). Here we consider the 0.45” gaussian

coronograph with 2 hour integrations at the MMT.

These simulations were used to define an optimal survey sample of 80 M0-F0 stars within 20 pc and less than 1 Gyr old. With this sample, and assuming $`dn/daa^{0.7}`$ and $`dn/dMM^{0.7}`$, we expect to detect 14 $`\pm `$ 3 companions with 5-15 $`M_{Jupiter}`$ at 13-50AU.

We do not necessarily expect our simple extrapolation to be valid given that at very small separations planets may have very different evolutionary histories than their larger separations counterparts (e.g. migration Trilling, Lunine, & Benz (2002)). Our survey will provide an indication of the shape of the separation and mass distribution functions based on the number of companions detected.

Figure 2 indicates our chances of detecting various numbers of planets given $`dn/daa^{0.7}`$ as well as three different mass distribution functions of the form $`dn/dMM^\alpha `$ with $`\alpha `$=\[-1,0,1\]. These three distribution functions are clearly distinguishable for the given sample. Assuming the given separation distribution function, we would be able to rule out $`dn/dMM^1`$ at 3$`\sigma `$ and $`dn/dMM^1`$ at 11$`\sigma `$ in the event of no planet detections.

The same process can be repeated for various separation distributions. Figure 6 indicates our chances of detecting various numbers of planets, while this time keeping $`dn/dMM^{0.7}`$ and varying the separation distributon function from $`dn/daa^\beta `$ $`\beta `$=\[-0.6,-0.2,0.1,0.7,2\]. For a null result we can rule out $`\beta `$=-0.2 at 2$`\sigma `$ and $`\beta `$=0.7 at 4$`\sigma `$.

These limits will provide valuable constraints on theories of planet formation and evolution as well as help focus the next generation of direct detection surveys.

We would like to thank IRLabs for their expert help, especially Elliott Solheid and Ken Salvestrini for designing the instrument dewar and testing the detector. Thanks also to Michael Lesser for helping us integrate AZCam with Klio. M.Freed acknowledges support from the NASA Graduate Student Researchers Program (NGT5-50394). We also acknowledge support from AFOSR as well as the NASA Astrobiology Institute “Laplace Center” node at the UofA.


text:
# Colloidal gelation, percolation and structural arrest

## Abstract

By means of molecular dynamics, we study a model system for colloidal suspensions where the interaction is based on a competition between attraction and repulsion. At low temperatures the relaxation time first increases as a power law as a function of the volume fraction and then, due to the finite lifetime of the bonded structures, it deviates from this critical behavior. We show that colloidal gelation at low temperatures and low volume fractions is crucially related to the formation of spanning long living cluster, as we suggested in the last few years. Besides agreeing with experimental findings in different colloidal systems, our results shed new light on the different role played by the formation of long living bonds and the crowding of the particles in colloidal structural arrest. PACS: 82.70.Dd, 64.60.Ak, 82.70.Gg

In attractive colloids a rich phenomenology is found in the temperature - volume fraction plane edinburgh -attr\_exp . At high temperature a hard sphere glass transition occurs at a volume fraction, $`\varphi 0.57`$. By decreasing the temperature, at high volume fraction, the effect of the short range attraction produces an attractive glass line with a reentrant line well described in the framework of the mode coupling theory attr ; attr\_sim , and confirmed by experiments. At low temperature and increasing the volume fraction a much more complex situation arises weitzclu , gron -sciortino3 , characterized by a cluster phase followed by a kinetic arrest with viscoelastic properties very similar to those found in polymer gelation.

One of the important problems is to understand the mechanism which gives rise to colloidal gelation. When phase separation does not occur, colloidal gelation has been interpreted in terms of cluster mode coupling theory kroy  or Wigner glass sciortino2 . In contrast to these approaches, we have suggested that at low volume fraction and low temperature, percolation phenomena play a crucial role in colloidal gelation characterized by the presence of a spanning cluster of bonded particles noi1 ; noi3  as also found in recent experiments bartlett . Of course, at high enough temperature and high volume fraction, one expects a structural arrest due essentially to the hard sphere part of the potential. At intermediate values of temperature and volume fractions, molecular dynamical studies have shown the occurrence of colloidal gelation phenomena characterized by a structural arrest, on average well described by the mode coupling theory of attractive glasses cates . Interestingly, the system is dynamically heterogeneous: The set of all particles can be partitioned in two subsets made of slow and fast particles and structural arrest occurs close to the onset of a spanning cluster of only the subset of slow particles cates . One may ask the question whether the same scenario for colloidal gelation can be smoothly extended from this intermediate regime to the low volume fraction region. In this letter, by using a molecular dynamics simulations of a DLVO-type potential dlvo  for charged colloidal systems, we give evidence of a clear crossover from the low volume fraction regime to the intermediate regime. In the first regime we find the bond lifetime much larger than the structural relaxation time. As a consequence, dynamically, the clusters behave as made of permanent bonds as in chemical gelation, and the regime is dominated by the dynamics of such clusters. Increasing $`\varphi `$, when these two time scales become comparable, the structural relaxation begins to be affected also by the crowding of the particles, and a clear crossover to a new glassy regime is found. Interestingly enough such crossover has also been found in some micellar systems at rather high volume fraction mallamaceetal ; noi1 , suggesting a common mechanism in both systems.

We consider a system of $`N=10000\varphi `$ particles, with $`\varphi =0.08`$$`0.25`$, interacting via a DLVO-type potential, which contains a Van der Waals type interaction plus an effective repulsion due to the presence of charges:

$$V(r)=ϵ\left[a_1\left(\frac{\sigma }{r}\right)^{36}a_2\left(\frac{\sigma }{r}\right)^6+a_3e^{\lambda (\frac{r}{\sigma }1)}\right],$$
(1)
where $`a_1=2.3`$, $`a_2=6`$, $`a_3=3.5`$, and $`\lambda =2.5`$ note\_yukawa . With these parameters the repulsion term dominates the Van der Waals attraction at long range, providing a short range attraction and a long range repulsive barrier. The effective repulsion in the potential prevents the liquid-gas phase separation and stabilize the size of the clusters, as expected gron . The potential is truncated and shifted at a distance of 3.5$`\sigma `$. To mimic the colloidal dynamics, we performed molecular dynamics simulations at constant temperature. Equations of motion were solved in the canonical ensemble (with a Nosé-Hoover thermostat) using a velocity Verlet algorithm nose  with a time step of $`0.001t_0`$ (where $`t_0=\sqrt{\frac{m\sigma ^2}{ϵ}}`$ and $`m`$ is the mass of the particles). We equilibrate the system at temperatures $`k_BT=0.2,0.23,0.25ϵ`$ and $`\varphi `$ increasing from $`0.07`$ to $`0.23`$ lamellar .

Here we analyze the structure and the dynamics by varying the temperature and volume fraction. We calculate first the static structure factor, $`S(k)`$, and the cluster size distribution, $`n(s)`$, where two particles are linked if their relative distance is smaller than the local maximum of the potential notahill .

At low temperatures, i.e. when the kinetic energy becomes smaller than the repulsive barrier, and low volume fractions, the static structure factor, $`S(k)`$, displays a peak around $`k_02`$ (see Fig. 1). This feature is due to the fact that the competing short range attraction and long range repulsion produce stable clusters of a typical size. The presence of a typical size is also manifested in the cluster size distribution, $`n(s)`$, plotted in Fig. 2 at $`k_BT=0.2ϵ`$ and $`\varphi =0.13`$. It is interesting to observe the secondary peak due to the fusion of two stable clusters. By increasing the volume fraction $`\varphi `$, a spanning cluster appears at $`\varphi =0.16`$ for $`k_BT=0.2ϵ`$, and the cluster size distribution typically displays a power law decay with an exponent $`\tau 2.18`$, consistent with the random percolation model staul . In Fig. 3 we show the dependence of the cluster size, $`s`$, on its radius, $`r`$, at $`\varphi _p`$: The data correspond to a compact structure for clusters of dimension $`s10`$ and radius $`r1`$, a fractal dimensionality $`d_f1.2`$ on intermediate length scales, and a crossover to $`d_f2.5`$ at larger length scales. These results suggest that at low volume fraction compact stable clusters form with typical size $`s10`$ and radius $`r1`$. By increasing the volume fraction a residual attractive interaction between the clusters produces tube-like structures with fractal dimension $`d_f1.2`$ up a size $`s60`$. By further increasing the volume fraction these structures coalesce to build a random percolating network. For the temperature here studied the percolation threshold $`\varphi _p`$ is weakly dependent on the temperature, $`0.16\varphi _p0.17`$.

A preliminary study of this model was published in Ref. noi3  and these results were presented in the meeting “Dynamical arrest”, hold in Bad Gastein in January 2005. This analysis has been lately developed in Ref. sciortino3 , where using a qualitatively similar interaction potential the authors have found again the behavior shown in Fig.s 2 and 3, and stressed the role of percolation. Interestingly they have also shown that this model is actually able to reproduce the local order typically observed in the experimental system bartlett .

In order to study the relaxation dynamics, we calculate different time autocorrelation functions at equilibrium. In Fig. 4 the incoherent scattering functions, $`F_s(k,t)`$, are plotted for $`k1.88`$, $`k_BT=0.25ϵ`$ and different $`\varphi `$. From these functions we extract a characteristic relaxation time, $`\tau _\alpha `$, defined by $`F_s(k,\tau _\alpha )=0.1`$, plotted in Fig. 6.

We also study the time correlation of the clusters in terms of the time autocorrelation functions of the bonds, $`B(t)`$, which are plotted in Fig. 5.

The function $`B(t)`$ is defined as

$$B(t)=\frac{_{ij}\left[n_{ij}(t)n_{ij}(0)n_{ij}^2\right]}{_{ij}\left[n_{ij}n_{ij}^2\right]},$$
(2)
where $`n_{ij}(t)=1`$ if particles $`i`$ and $`j`$ are linked at time $`t`$, $`n_{ij}(t)=0`$ otherwise. The data show that, differently from the behavior of $`F_s(k,t)`$, the bond relaxation does not change much below the percolation threshold, $`\varphi _p`$, and instead it becomes slower and slower above the percolation threshold. By fitting the curves with stretched exponential functions, we obtain the bond lifetime $`\tau _b`$. In Fig. 6 $`\tau _b`$ is compared with $`\tau _\alpha `$. The data show that, at low volume fractions, the bond lifetime is much larger than the density relaxation time and in this region $`\tau _\alpha `$ can be fitted by a power law. This result can be understood by considering that the system dynamically behaves as made of permanent clusters and the structural relaxation is directly related to the size of the clusters, like in chemical gelation. At higher volume fraction close to the percolation threshold, the relaxation time and the bond lifetime eventually become of the same order of magnitude and one observes a deviation from the power law critical behavior.

Within this scenario, the percolation of long-living structures is the mechanism which produces the enhancement of the structural relaxation time at low volume fraction: Close to the percolation threshold, $`\tau _\alpha `$ increases fast enough to reach the scale of the bond lifetime, and the dynamics crosses over to the glassy regime.

In fact in this region on the time scale of the structural relaxation time, the bonds cannot be considered as permanent anymore. This explains the crossover to a glassy regime, in which the structural relaxation time depends on the bond finite lifetime and the crowding of the particles.

At higher temperature when the bond lifetime does not play any role, structural arrest is only due to the hard sphere component of the potential.

By decreasing the temperature we find that the critical regime becomes more extended, due to the increase of the bond lifetime (see Fig. 7), and the critical volume fraction $`\varphi _c`$ tends to the percolation threshold $`\varphi _p`$. Therefore at very low temperatures we expect the structural arrest to coincide with the onset of a spanning cluster on the time scale of the numerical simulations, in agreement with previous results noi3  and recent experimentsbartlett .

In order to further understand and characterize the crossover region, we have calculated the dynamic susceptibility, defined as the fluctuations of $`F_s(k,t)`$ toninelli :

$$\chi _4(k,t)=N(F_s(k,t)^2F_s(k,t)^2).$$
(3)
Similar quantities, largely studied in glassy systems (see Refs. franz ; glotzer ), are related to the presence of dynamic heterogeneities, spatial correlations between particles that move faster than the average. In Fig. 8, $`\chi _4(k,t)`$ is plotted as function of $`t`$, for $`k1.88`$, $`k_BT=0.25ϵ`$ and different values of $`\varphi `$. We find that $`\chi _4(k,t)`$ exhibits a maximum at time $`t^{}`$, which first increases and then decreases signalling the presence of a crossover as clearly shown in Fig. 9, where the time, $`t^{}`$, and the maximum value, $`\chi _4(k,t^{})`$, are plotted as a function of the volume fraction: At low volume fraction they both follow a power law diverging, whereas a deviation from this critical behavior appears at volume fraction higher than the percolation threshold.

The presence of a maximum at a time $`t^{}`$, increasing together with the height of the maximum was recently found in experiments on colloidal gels as colloidal gelation was approached cipelletti .

In conclusions in our model the presence of two well separated time scales at low temperature and low volume fraction originates a colloidal gelation related to the formation of persistent structures: Increasing the volume fraction the two time scales become comparable and a crossover from a gel-like behavior to a glass-like one is found.

This work has been partially supported by the Marie Curie Fellowship HPMF-CI2002-01945 and Reintegration Grant MERG-CT-2004-012867, EU Network Number MRTN-CT-2003-504712, MIUR-PRIN 2004, MIUR-FIRB 2001, CRdC-AMRA, INFM-PCI.