text:
# Fireballs, Flares and Flickering

## 1. Introduction

AE Aqr is a long period CV ($`P_{\mathrm{orb}}=9.88\mathrm{h}`$), implying that the secondary is somewhat evolved. Coherent oscillations in optical (Patterson 1979), ultraviolet (Eracleous et al. 1994), and X-ray (Patterson et al. 1980) lightcurves reveal the white dwarf’s 33s spin period. The oscillation has two unequal peaks per spin cycle, consistent with broad hotspots on opposite sides of the white dwarf $`15^{}`$ above and below the equator (Eracleous et al. 1994). These oscillations are strongest in the ultraviolet, where their spectra show a blue continuum with broad Ly $`\alpha `$ absorption consistent with a white dwarf ($`\mathrm{log}g8`$) atmosphere with $`T3\times 10^4\mathrm{K}`$ (Eracleous et al. 1994). The simplest interpretation is accretion heating near the poles of a magnetic dipole field tipped almost perpendicular to the rotation axis.

An 11-year study of the optical oscillation period (de Jager et al. 1994) showed that the white dwarf is spinning down at an alarming rate. Something extracts rotational energy from the white dwarf at a rate $`I\omega \dot{\omega }60\nu L_\nu `$ ie. some 60 times the luminosity of the system. We now believe that to be a magnetic propeller. This model was first proposed by Wynn, King, & Horne (1995) and Eracleous & Horne (1996) and expanded upon in Wynn, King & Horne (1997) with comparison of observed and modelled tomograms. Further work on the flaring region was reported in Welsh, Horne & Gomer (1998) and Horne (1999).

The gas stream emerging from the companion star through the L1 nozzle encounters a rapidly spinning magnetosphere. The rapid spin makes the effective ram pressure so high that only a low-density fringe of material becomes threaded onto field lines. Most of the stream material remains diamagnetic, and is dragged toward co-rotation with the magnetosphere. As this occurs outside the co-rotation radius, this magnetic drag propels material forward, boosting its velocity up to and beyond escape velocity. The material emerges from the magnetosphere and sails out of the binary system. This process efficiently extracts energy and angular momentum from the white dwarf, transferring it via the long-range magnetic field to the stream material, which is expelled from the system. The ejected outflow consists of a broad equatorial fan of material launched over a range of azimuths on the side away from the secondary..

The material stripped from the gas stream and threaded by the field lines has a different fate, one which we believe gives rise to the radio and X-ray emission. This material co-rotates with the magnetosphere while accelerating along field lines either toward or away from the white dwarf under the influences of gravity and centrifugal forces. The small fraction of the total mass transfer that leaks below the co-rotation radius at $`5R_{\mathrm{wd}}`$ accretes down field lines producing the surface hotspots responsible for the 33s oscillations. Particles outside co-rotation remain trapped long enough to accelerate up to relativistic energies through magnetic pumping, eventually reaching a sufficient energy density to break away from the magnetosphere (Kuijpers et al. 1997). The resulting ejection of balls of relativistic magnetized plasma is thought to give rise to the flaring radio emission (Bastian, Dulk & Chanmugam 1998a,b).

In many studies the lightcurves exhibit dramatic flares, with 1-10 minute rise and fall times (Patterson 1979; van Paradijs, Kraakman & van Amerongen 1989; Bruch 1991; Welsh, Horne & Oke 1993). The flares seem to come in clusters or avalanches of many super-imposed individual flares separated by quiet intervals of gradually declining line and continuum emission (Eracleous & Horne 1996; Patterson 1979). These quiet and flaring states typically last a few hours. Power spectra computed from the lightcurves have a power-law form, with larger amplitudes on longer timescales. Similar red noise power spectra are seen in active galaxies, X-ray binaries, and other cataclysmic variables, and is therefore regarded as characteristic of accreting sources in general. However, flickering in other cataclysmic variables typically has an amplitude of 5-20% (Bruch 1992), contrasting with factors of several in AE Aqr. If the mechanism is the same, then it must be weaker or dramatically diluted in other systems. The extreme behaviour in AE Aqr may thus allow us to probe the underlying physics involved where the resultant effects are clearest.

The optical and ultraviolet spectra of the AE Aqr flares are not understood at present except in the most general terms. The lines and continua rise and fall together, with little change in the equivalent widths or ratios of the emission lines (Eracleous & Horne 1996). This suggests that the flares represent changes in the amount of material involved more than changes in physical conditions. Ultraviolet spectra from HST reveal a wide range of lines representing a diverse mix of ionization states and densities. Eracleous & Horne (1996) concluded that a large fraction of the line-emitting gas the density is in the range $`n10^{15}`$$`10^{17}\mathrm{m}^3`$ and that denser regions likely also exist. The CIV emission is unusually weak, suggesting non-solar abundances consistent with the idea that the secondary is evolved. Latest results of ultraviolet observations and abundance analysis of the secondary are presented by Mouchet (this volume) and Bonnet-Bidaud (this volume).

What mechanism triggers these dramatic optical and ultraviolet flares? Clues come from multi-wavelength co-variability and orbital kinematics. Simultaneous VLA and optical observations show that the radio flux variations occur on similar timescales but are not correlated with the optical and ultraviolet flares, which therefore require a different mechanism (Abada-Simon et al. 1995). It has been proposed that the flares represent modulations of the accretion rate onto the white dwarf, so that they should be correlated with X-ray variability. Some correlation was found, but the correlation is not high. However, HST observations discard this model, because the ultraviolet oscillation amplitude is unmoved by transitions between the quiet and flaring states (Eracleous & Horne 1996). This disconnects the origins of the oscillations and flares, and the oscillations arise from accretion onto the white dwarf, so the flares must arise elsewhere.

Further clues come from emission line kinematics. The emission line profiles may be roughly described as broad Gaussians with widths $`1000`$ km s-1, though they often exhibit kinks and sometimes multiple peaks. Detailed study of the Balmer lines (Welsh et al. 1998) indicates that the new light appearing during a flare can have emission lines shifted from the line centroid and somewhat narrower, $`300`$ km s-1. Individual flares therefore occupy only a subset of the entire emission-line region.

The emission line centroid velocities vary sinusoidally with orbital phase, with semi-amplitudes $`200`$ km s-1 and maximum redshift near phase $``$0.8 (Welsh et al. 1998). These unusual orbital kinematics are shared by both ultraviolet and optical emission lines. The implication is that the flares arise from gas moving with a $`200`$ km s-1 velocity vector that rotates with the binary and points away the observer at phase 0.8. This is hard to understand in the standard model of a cataclysmic variable star, though many (eg. SW Sex systems) show similar anomalous emission-line kinematics (Thorstensen et al. 1991) suggesting a common underlying mechanism. In particular, Kepler velocities in an AE Aqr accretion disc would be $`>600`$ km s-1 (though we believe no disc to be present.) The gas stream has a similar direction but its velocity is $`1000`$ km s-1. A success of the magnetic propeller model is its ability to account for the anomalous emission-line kinematics. The correct velocity amplitude and direction occurs in the exit fan just outside the Roche lobe of the white dwarf. But the question remains of why the flares are ignited here, several hours after the gas slips silently through the magnetosphere.

The key insight which solved this puzzle was the realization that the magnetic propeller acts as a blob sorter. More compact, denser, diamagnetic blobs are less affected by magnetic drag. They punch deeper into the magnetosphere and emerge at a larger azimuth with a smaller terminal velocity. These compact blobs can therefore be overtaken by ‘fluffier’ blobs ejected with a larger terminal velocity in the same direction, having left the companion star somewhat later but having spent less time in the magnetosphere (Wynn et al. 1997). The result is a collision between two gas blobs, which can give rise to shocks and flares. Calculations of the trajectories of magnetically propelled diamagnetic blobs with different drag coefficients indicate that they cross in an arc-shaped region of the exit stream, in just the right place to account for the orbital kinematics of the emission lines (Welsh et al. 1998; Horne 1999). Figure 1 shows how these trajectories map to a locus of points in the lower left quadrant of the Doppler map that is otherwise difficult to populate.

There remains the problem of developing a quantitative understanding of the unusual emission-line spectra resulting from the collision.

## 2. Fireball Models

### 2.1. Typical Parameters

We can derive rough estimates for the physical parameters associated with the colliding blobs by using the typical rise time $`t300\mathrm{s}`$ and an observed optical flux $`f_\nu 50`$ mJy with a closing velocity $`V300\mathrm{km}\mathrm{s}^1`$ for the blobs, together with a distance of 100pc. A summary of typical values is provided in Table 1.

### 2.2. Expansion

We envision an expanding fireball emerging from the aftermath of a collision between two masses $`m_1+m_2=M`$. With nothing to hold it back, the hot ball of gas expands at the initial sound speed, launching a fireball. We adopt a uniform, spherically symmetric, Hubble-like expansion $`V=Hr_0=Hr/\beta `$, in which the Eulerian radial coordinate $`r`$ of a gas element is given in terms of its initial position $`r_0`$ and time $`t`$ by

$$r(r_0,t)=r_0+v(r_0)t=r_0+Hr_0t=r_0(1+Ht)r_0\beta .$$
(1)
This defines an expansion factor $`\beta `$ which we can use as a dimensionless time parameter. The ‘Hubble’ constant is set by the initial conditions $`Ha_\mathrm{i}v(a_\mathrm{i})c_{\mathrm{s},\mathrm{i}}`$. If we can determine parameters at some time $`t=0`$ ie. $`\beta =1`$ for the lengthscale ($`a_0`$), temperature ($`T_0`$) and mass ($`M`$), we can derive the time evolution as follows. We adopt a Gaussian density profile

$$\rho (r,t)=\rho _0\beta ^3e^{\eta ^2}$$
(2)
where $`\rho _0`$ is the central density at $`t=0`$ and $`\eta \frac{r_0}{a_0}=\frac{r}{a}`$ is a dimensionless radius coordinate: the radius $`r`$ scaled to the lengthscale $`a=\beta a_0`$. This Gaussian density is motivated by the Gaussian shapes of observed, Doppler-broadened, line profiles.

### 2.3. Cooling

We consider two cooling schemes: adiabatic and isothermal. The adiabatic fireball cools purely as a result of its expansion and corresponds to a situation where the radiative and recombination cooling rates are negligible. In contrast, the isothermal model maintains a fixed temperature throughout its evolution. A truly isothermal fireball would require a finely balanced energy source to counteract the expansion cooling. However, it may be an appropriate approximation to a situation where a photospheric region dominates the emission and presents a fixed effective temperature to the observer as a result of the stong dependence of opacity on temperature. For adiabatic expansion we can show that an initially uniform temperature distribution remains uniform and so,

$$T(r,t)=T_0\beta ^{3(\gamma 1)}.$$
(3)
With $`\gamma =5/3`$ for a monatomic gas, $`T\beta ^2`$, and so the sound speed $`c_\mathrm{s}T^{1/2}\beta ^1`$. The sound crossing time $`a/c_\mathrm{s}\beta ^2`$. and so the fireball becomes almost immediately supersonic.

### 2.4. Ionization Structure

In the LTE approximation, the density and temperature determine the ionization state of the gas at each point in space and time through the solution of a network of Saha equations. Atomic level populations are similarly determined through Boltzmann factors and partition functions.

Once we have determined the evolution of $`T`$ and $`\rho `$ with time, we can follow the evolution of the ionization structure. At a given time, the various ionic species form onion layers of increasing ionization on moving outwards from the fireball centre. The adiabatic temperature evolution causes the ionization states to alter rapidly throughout the structures when the fireball passes the appropriate critical temperatures. Isothermal evolutions show the gentler dependance of ionization on density. The spatial structure remains roughly constant and only evolves slowly as the density drops.

## 3. Radiative Transfer

The radiative transfer equation has a formal solution

$$I=Se^\tau 𝑑\tau ,$$
(4)
where $`I`$ is the intensity of the emerging radiation, $`S`$ is the source function, and $`\tau `$ is the optical depth measured along the line of sight from the observer. This integral sums contributions $`Sd\tau `$ to the radiation intensity, attenuating each by the factor $`e^\tau `$ because it has to pass through optical depth $`\tau `$ to reach the observer.

Since we have assumed LTE, the source function is the Planck function, and opacities both for lines and continuum are also known once the velocity, temperature, and density profiles and element abundances are specified. The integral can therefore be evaluated numerically, either in this form or more quickly by using Sobolev resonant surface approximations.

The above line integral gives the intensity $`I(y)`$ for lines of sight with different impact parameters $`y`$. We let $`x`$ measure distance from the fireball centre toward the observer, and $`y`$ the distance perpendicular to the line of sight. The fireball flux, obtained by summing intensities weighted by the solid angles of annuli on the sky, is then

$$f(\lambda )=_0^{\mathrm{}}I(\lambda ,y)\frac{2\pi y}{d^2}𝑑y$$
(5)
where $`d`$ is the source distance.

## 4. Comparison with Observations

### 4.1. Observed Flare in AE Aqr

To test the fireball models, we compared them to high time-resolution optical spectra of an AE Aqr flare taken with the Keck telescope (Skidmore et al. 2003). The lightcurve formed from these data is plotted in Fig. 2. Although the dataset only lasts around 13 minutes it shows the end of one flare and beginning of another. The major features of the flaring observed in the system are apparent in the lightcurve.

We extracted the spectrum of the second flare by subtracting the quiescent spectrum (Q) from the spectrum at the top of the new flare (P). The resultant optical spectrum is plotted in the bottom panel of Fig. 3.

### 4.2. Spectral Fit at Peak of Flare

To estimate the fireball parameters $`M`$, $`T`$ and $`a`$ from the observed optical spectra we began with consideration of the continuum flux. We considered both Population I and Population II abundances. Population I abundances were set using the solar abundances presented by (Däppen 2000). For population II composition, we decreased from 0.085 to 0.028 the ratio of the helium to hydrogen number densities and reduced the ratio for each metal abundance to 5% of solar (Bowers & Deeming 1984). In both cases we used an amoeba algorithm (Press et al. 1986) to optimize the fit by reducing $`\chi ^2`$ to a minimum for 6 optical continuum fluxes. The best fit parameters $`M`$, $`T`$ and $`a`$ and related secondary quantities are listed in Table 2.

Both compositions arrive at a model fit with remarkably similar calculated fluxes which have a slightly steeper slope to the Paschen continuum than the observations. The fits also produce a large Balmer jump but appears to have difficulty in creating one quite as strong as that observed. The central densities implied by both sets of parameters are consistent with the typical values for the gas stream in section 2.1. and with the temperature range and emitting area of Beskrovnaya et al. (1996).

We derived the expansion rate from consideration of the H$`\beta `$ and nearby HeI equivalent widths. In a fireball model, increasing the expansion velocity spreads out the line opacity and thereby increases the equivalent widths of the emission lines as long as they remain optically thick. The H$`\beta `$ region was chosen because the observed H$`\alpha `$ line shows a complex structure. We blurred the model spectra with a Gaussian of $`9.8\text{Å}`$ FWHM and assumed a distance of $`102`$ pc (Friedjung 1997). The H$`\beta `$ and HeI line equivalent widths suggest expansion velocities of $`240\text{km}\text{s}^1`$ and $`100\text{km}\text{s}^1`$ respectively at $`a_0`$. Consequently, we adopted an expansion velocity at $`a_0`$ of $`170\text{km}\text{s}^1`$ for the simulations.

The optical spectra produced with both compositions are compared to the observed mean spectrum at the flare peak in Fig. 3. We see that the parameters, derived from the peak continuum fluxes alone, allow us to derive optical spectra with integrated line fluxes and ratios comparable to the observations. Both models have line widths narrower than the observations, and hence peak line fluxes greater than those observed, which may result from additional unidentified instrumental blurring. The Population I models exhibit a flat Balmer decrement of saturated Balmer lines that is also apparent in the data. Comparing the lines that are present, however, particularly in the $`40005000`$ Å range, the Population II model appears to produce a much better fit: suppressing the metal lines which do not appear in the observed spectra.

### 4.3. Lightcurves and Spectral Evolution

The timescale for the evolution of the fireball implied by the observed lightcurve can be reproduced by the parameters derived solely from fitting to the peak spectrum. The evolutionary timescale $`a_0/v(a_0)=1/H`$ is approximately $`300`$ s and $`560`$ s for Population I and II parameters respectively (Table 2). The encouraging result that these values are similar to the observed timescale suggests that our fireball models are on the right track.

Fig. 4 shows the lightcurves and spectral evolution for Population II abundances derived for isothermal and adiabatic evolutions (‘radiative’ models will be discussed later). In each case, the model fireball is constrained to evolve through the state fitted to the observations at the peak of the flare.

The isothermal models produces a lightcurve that rises and falls, as observed. The isothermal fireball spectra, both on the rise and on the fall of the flare, exhibit strong Balmer emission lines and a Balmer jump in emission. Thus the isothermal fireball model provides a plausible fit not only to the peak spectrum but also to the time evolution of the AE Aqr flare.

The adiabatic model fails to reproduce the observed lightcurve or spectra. The adiabatic model matches the observed spectrum at one time, but is too hot at early times and too cool at late times. As a result the lightcurve declines monotonically from early times rather than rising to a peak and then falling, and the spectra have lines that do not match the observed spectra. This disappointing performance is perhaps not too surprising, since the temperature in this model falls through a wide range while the observations suggest that the temperature is always around $`10^4`$ K. To explain the rise phase of the observed lightcurve, for an adiabatic evolution, we must interpret it as resulting from the collision of the initial gas blobs.

## 5. Discussion

The spectra produced with Population II composition reproduce the optical observations well. The fitted parameters are consistent with both the expected conditions in the mass transfer stream and with the lower limits on the density implied by the ultraviolet observations ($`n10^{15}\text{m}^3`$). The presence of high ionization and semi-forbidden lines in the ultraviolet data may be understood qualitatively in terms of their formation in the low density outer fringes of the fireball and the wide variety of ionization states by the wide range of densities in the expanding fireball structure. For $`\beta =1`$, mean molcular weight $`\mu =0.53`$ and Population II parameters, the density of $`10^{17}\text{m}^3`$, important for the ultraviolet semi-forbidden lines, occurs at $`\eta 2.3`$. We can show that this lies outside the limit for LTE behaviour (Pearson, Horne & Skidmore 2003). Quantitative fits for the fringe and ultraviolet behaviour, therefore, remain to be addressed in the non-LTE regime.

The behaviour of the lightcurve for an isothermal model is much more similar to the observed lightcurve than the adiabatic model, giving a peak close to the observed time without the need to invoke the collision process to create a rising phase. Clearly though, an expanding gas ball would be expected to cool both from radiative and adiabatic expansion effects.

In short, we find that isothermal fireballs reproduce the observations rather well, whereas adiabatically cooling fireballs fail miserably! How can the expanding gas ball present a nearly constant temperature to the observer?

### 5.1. Thermostat

We examine two possible mechanisms which may operate to maintain the apparent fireball temperature in the $`1`$$`2\times 10^4\text{K}`$ region. The first relies on the fact that both free-free and bound-free opacity decrease with temperature and thus opacity peaks just above the temperature at which hydrogen begins to recombine. In this model a hot core region is surrounded by a cooler blanket of material. Radiating at the hotter temperature, the core would rapidly lose its thermal energy through radiative cooling. However, we can envision a situation where the material just outside the core has reached the high opacity regime and is absorbing significant amounts of energy from the core’s radiation field. This would help to counteract the effect of adiabatic cooling in this shell and may maintain an effectively isothermal blanket around the core in the following way. If the temperature in the blanket were to rise, the opacity would decrease, more energy would escape and the blanket would cool again. Similarly, until recombination, cooling would result in higher opacity and therefore greater heating of the blanket region. The core region would thus only be cooling through adiabatic expansion. For such a thermostatic mechanism to work, the optical depth through the blanket region would need to be high which would be consistent with it being the photosphere as seen by an outside observer.

Simple models have been attempted to mimic the effect of a photosphere shrinking due to radiative cooling but these have been unsuccesful in reproducing the observations (Pearson et al. 2003). This occurs because the cooling front occurs well outside the photosphere. The thin shell of cooling material at around $`10^4\text{K}`$ does not provide sufficient opacity to mimic the isothermal behaviour and we therefore “see” into the adiabatically-cooling core. The radiative model thus shows very similar behaviour to the adiabatic models.

Realistic simulations of such a model require a time-dependent radiation hydrodynamical treatment of the radiation field and its heating effect; a more sophisticated method than the one so far employed.

### 5.2. External Photoionization

Alternatively, the ionization of the fireball might be held at a temperature typical of the $`1`$$`2\times 10^4\text{K}`$ range through photoionization by the white dwarf. In Fig. 5, we plot the hydrogen ionization and recombination rates as a function of temperature for a typical electron number density $`n_\mathrm{e}=3\times 10^{18}\mathrm{m}^3`$ using routines for photoionization from Verner et al. (1996), radiative recombination (Verner & Ferland 1996), three-body recombination (Cota 1987) and collisional ionization by Verner using the results of Voronov (1997). We calculate a conservative overestimate for the self-photoionization of the fireball assuming the sky is half filled by blackbody radiation at the given fireball temperature. Also shown is the photoionization due to a blackbody with the same radius and distance as the white dwarf. The general white dwarf temperature is uncertain but, as mentioned earlier, there is evidence to suggest that the hot spots have a temperature of around $`3\times 10^4\text{K}`$. We can see that the flux from the white dwarf overtakes that from the fireball itself as the dominant ionization mechanism at $`10^4\text{K}`$. We have assumed here a typical white dwarf to fireball distance equal to that of the L1 point from the white dwarf (although in a different direction). The recombination rates for this plot have been calculated assuming $`n_\mathrm{e}`$ is a constant and so for temperatures below about $`10^4\text{K}`$ are conservative overestimates of the LTE recombination rate that will occur in our fireball. In spite of this, the white dwarf photoionization at $`3\times 10^4\text{K}`$ comfortably exceeds the recombination rates down to at least $`10^3\text{K}`$ and, hence, we would expect the fireball to remain almost completely ionized in this case. This suggests that the white dwarf radiation field may be the cause of the fireball appearing isothermal when we consider the lightcurve behaviour and consistent line strengths and ratios. Ionization by an external black body at a different temperature will clearly cause our fireball to no longer be in LTE and require a non-LTE model to follow the fireball behaviour accurately.

## 6. Summary

We have shown for AE Aqr how the observed flare spectrum and evolution is reproducible with an isothermal fireball at $`T18000`$ K with Population II abundances but not when adiabatic cooling is incorporated. We suspect that the cause of the apparently isothermal nature is a combination of two mechanisms. First, a nearly isothermal photosphere which is self-regulated by the temperature dependance of the continuum opacity and the hydrogen recombination front and second, particularly at late times, by the ionizing effect of the white dwarf radiation field.

The observed photosphere in the isothermal models is initially advected outward with the flow before the decreasing density causes the opacity to drop and the photosphere to shrink to zero size. This gives rise to a lightcurve that rises and then falls. Emission lines and edges arise because the photospheric radius is larger at wavelengths with higher opacity. Surfaces of constant Doppler shift are perpendicular to the line of sight and thus Gaussian density profiles give rise to Gaussian velocity profiles.

In a purely LTE fireball we expect high ionization states to occur in the outer regions as a result of the lower density. In reality, of course, once the LTE boundary is crossed, the ionization will drop as the net recombination rate gradually ekes away at the ions and we make a transition towards coronal equilibrium. In principle, we can follow the non-LTE evolution of the gas through this outer region and beyond the ion-electron equilibrium radius to the non-equilibrium ionization states present in the outer fringe. The higher ionization would then result from the higher temperature of the fireball when the region in question crossed the LTE boundary. In practice, we expect that this correction may have little effect on the optical spectra because the emission is dominated by the higher density inner regions of the fireball. This effect may well become important, however, when we consider the ultraviolet spectra where high ionization lines are present; and which are not reproduced by the current models.

Improved modelling of these fireballs offers us the chance to probe the chemical composition of the secondary star in AE Aqr in a complimentary regime to those presented by Bonnet-Bidaud (this volume). Preliminary results from SS Cyg encourage us that the fireball models will be applicable more generally to the flickering observed in most CV systems. In such ‘disc flickering’ systems, the fireballs may be launched by local heating from magnetic reconnection events rather than from blob-blob collisions. Similar spectral features are present in the flickering data (see Skidmore et al. in this volume) with similar temporal evolution to the AE Aqr flares. Initial fitting using the same procedure as for the AE Aqr observations give good results that suggest that a common mechanism may indeed be at work.

### Acknowledgements

KJP would like to thank PPARC, NSF grant 115-30-5177 and the IAU whose financial support made attendance at the colloquium possible.

## References

Abada-Simon, M., Bastian, T. S., Horne, K. D., Robinson, E. L., Bookbinder, J. A. 1995, in Buckley, D. A. H., Warner, B., eds., Proc. Cape Workshop on Magnetic Cataclysmic Variables, ASP Conf. Series, 355

Bastian, T. S., Dulk, G. A., Chanmugam, G. 1988a, ApJ, 324, 431

Bastian, T. S., Dulk, G. A., Chanmugam, G. 1988b, ApJ, 330, 518

Beskrovnaya, N., Ikhsanov, N., Bruch, A., Shakhovskoy, N. 1996, A&A, 307, 840

Bowers, R. L., Deeming, T. 1984, Astrophysics I Stars, Jones and Bartlett, Boston

Bruch, A. 1991, A&A, 251, 59

Bruch, A. 1992, A&A, 266, 237

Cota, S. A. 1987, PhD thesis, Ohio State Univ.

Däppen, W. 2000, in Cox, A. N., Allen’s Astrophysical Quantities, Chapter 3, Springer-Verlag, New York

Eracleous, M., Horne, K. D., Robinson, E. L., Zhang, E.-H., Marsh, T. R., Wood, J. H. 1994, ApJ, 433, 313

Eracleous, M., Horne, K. D. 1996, ApJ, 471, 427

Friedjung, M. 1997, NewA, 2, 319

Horne, K. D. 1999, in Hellier, K., Mukai, K., eds, Annapolis Workshop on Magnetic Cataclymic Variables, ASP Conf. Series, 157, 357

de Jager, O. C., Meintjes, P. J., O’Donoghue, D., Robinson, E. L. 1994, MNRAS, 267, 577

Kingdon, J. B., Ferland, G. H. 1996, ApJS, 106, 205

Kuijpers, J., Fletcher, L., Abada-Simon, M., Horne, K. D., Raadu, M. A., Ramsay, G., Steeghs, D. 1997, A&A, 322, 242

van Paradijs, J., Kraakman, H., van Amerongen, S. 1989, A&AS, 79, 205

Patterson, J. 1979, ApJ, 234, 978

Patterson, J., Branch, D., Chincarini, G., Robinson, E. L. 1980, ApJ, 240, L133

Pearson, K. J., Horne, K. D., Skidmore, W. 2003, MNRAS, 338, 1067

Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P. 1986, Numerical Recipes in Fortran, Cambridge Univ. Press

Skidmore, W., O’Brien, K., Horne, K. D., Gomer, R., Oke, J. B., Pearson, K. J. 2003, MNRAS, 338, 1057

Thorstensen, J. R., Ringwald, F. A., Wada, R. A., Schmidt, G. A., Norsworthy, J. E. 1991, AJ, 102, 272

Verner, D. A., Ferland, G. J. 1996, ApJS, 103, 467

Verner, D. A., Ferland, G. J., Korista, K. T., Yakovlev, D. G. 1996, ApJ, 465, 487

Voronov, G. S. 1997, Atomic Data and Nuclear Data Tables, 65, 1

Welsh, W. F., Horne, K., Oke, B. 1993, ApJ, 406, 229

Welsh, W. F., Horne, K. D., Gomer, R. 1998, MNRAS, 298, 285

Wynn, G. A., King, A. R., Horne, K. D. 1995, in Buckley, D. A. H., Warner, B., eds., Cape Workshop on Magnetic Cataclysmic Variables, ASP Conf. Series, 85, 196, Astron. Soc. Pacific., San Francisco

Wynn, G. A., King, A. R., Horne, K. D. 1997, MNRAS, 286, 436


format: text

text:
# Radial Velocity Studies of Close Binary Stars. VIII1footnote 11footnote 1Based on the data obtained at the David Dunlap Observatory, University of Toronto.

## 1 INTRODUCTION

This paper is a continuation in a series of papers of radial-velocity studies of close binary stars (Lu & Rucinski, 1999; Rucinski & Lu, 1999; Rucinski et al., 2000; Lu et al., 2001; Rucinski et al., 2001, 2002a) and presents data for the seventh group of ten close binary stars observed at the David Dunlap Observatory. Selection of the targets is quasi-random: At a given time, we observe a few dozen close binary systems with periods shorter than one day, brighter than 11 magnitude and with declinations $`>20^{}`$; we publish the results in groups of ten systems as soon as reasonable orbital elements are obtained from measurements evenly distributed in orbital phases. For technical details and conventions, and for preliminary estimates of errors and uncertainties, see the interim summary paper Rucinski (2002b, Paper VII).

This paper is structured in the same way as the previous papers in that most of the data for the observed binaries are in two tables consisting of the radial-velocity measurements (Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.) and their sine-curve solutions (Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.). Section 2 of the paper contains brief summaries of previous studies for individual systems and comments on the new data. Figures Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto. show the radial velocity data and solutions. While in the previous papers we showed only some most interesting broadening functions (BF), Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto. in this paper shows the BF’s for all systems; the functions have been selected from among the best defined ones around the orbital phase of 0.25.

The observations reported in this paper have been collected between March 2000 and April 2002; the ranges of dates for individual systems can be found in Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.. All systems discussed in this paper, except V523 Cas, have been observed for radial-velocity variations for the first time. We have derived the radial velocities in the same way as described in previous papers; see Paper VII for a discussion of the broadening-function approach used in the derivation of the radial-velocity orbit parameters: the amplitudes, $`K_i`$, the center-of-mass velocity, $`V_0`$, and the time-of-eclipse epoch, $`T_0`$.

The data in Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto. are organized in the same manner as in previous papers. In addition to the parameters of spectroscopic orbits, the table provides information about the relation between the spectroscopically observed epoch of the primary-eclipse T0 and the recent photometric determinations in the form of the $`OC`$ deviations for the number of elapsed periods E. It also contains our new spectral classifications of the program objects.

## 2 RESULTS FOR INDIVIDUAL SYSTEMS

### 2.1 V410 Aur

V410 Aur was discovered by the Hipparcos satellite. The somewhat sparsely covered Hipparcos light curve is rather typical for a W UMa-type system; it has an amplitude of about 0.33 mag. Total eclipses are entirely possible, given the large semi-amplitudes $`K_i`$ of our solution, suggesting an orbital inclination $`i90`$ degrees. The system is a spectroscopic triple, with the third component having the relative brightness of $`L_3/(L_1+L_2)=0.26\pm 0.01`$. The photometric amplitude corrected for the third light is 0.43 mag, which actually is large and inconsistent with the the mass-ratio $`q_{sp}=0.14`$; the maximum amplitude for the inclination of $`i=90`$ degrees is expected to be around 0.33 – 0.38 mag. (Rucinski, 2001).

The individual velocities of the third component did not show any trace of a “cross-talk” (i.e. did not correlate with the binary phase) although the third-star peak in the broadening function was strong and always projected into the broadening function signature of the primary component (see Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.) making the orbital solution somewhat poorer than for similar systems without companions (see Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.). The radial velocity observations of the third component are given in Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.. Our observations were made in two groups, around HJD 2,451,870 and HJD 2,452,280. The average velocities of the third component for these groups were $`V_3=48.38\pm 0.30`$ km s-1 and $`43.90\pm 0.43`$ km s-1; both differ significantly from the mean velocity of the close binary, $`V_0=36.94\pm 3.10`$ km s-1. Unfortunately, we do not have enough observations and they are not accurate enough to check if the difference in the average $`V_3`$ is reflected in a change of the binary center-of-mass velocity $`V_0`$; our solution for the binary uses one constant value for $`V_0`$.

The orbital ephemeris based on the Hipparcos data shows a large $`OC`$ deviation of 0.10 days. We arbitrarily reduced the $`OC`$ by shifting the time of the primary eclipse by that amount under an assumption that V410 Aur is a contact binary of the A type (i.e. that the hotter component is the more massive one).

The third body in the system was probably the reason for a large error in the Hipparcos parallax, $`4.77\pm 5.39`$ mas, so that we have no distance estimate for the system. There is a discrepancy between the $`V_{max}`$ estimated from the Hipparcos $`H_P=9.98`$ and from the Tycho–2 (Hog et al., 2000) photometry, $`V_{max}=10.18`$ (transformed from $`V_T`$). Our estimate of the spectral type, G0/2V, is consistent with the average $`(BV)=0.56`$ from Tycho–2.

### 2.2 V523 Cas

V523 Cas is the only system in this group which was observed spectroscopically before our observations (Milone et al., 1985). The system is an interesting one because – with its very short orbital period of 0.237 days, right next to CC Com at 0.221 days – it is very close to the still-unexplained abrupt cutoff in the orbital period distribution of contact binaries. Several photometric studies of V523 Cas have been recently rediscussed by Lister et al. (2000) where all previous references are given. The photometric solutions encountered the common problem of a featureless light curve poorly constraining the geometrical parameters, with the crucial parameter of the mass ratio taking a large range of values, with the photometric ones usually above the value determined by Milone et al. of $`q_{sp}=0.42\pm 0.02`$.

Our radial velocity results are well defined (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.). The broadening function (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.) very clearly shows the more compact and slightly hotter (or having a higher surface brightness) less-massive component. The disparity of the component shapes is quite striking in the broadening function; it seems as if the secondary were a detached component (but still quite large) in a semi-detached binary. Our spectroscopic mass ratio, $`q_{sp}=0.516\pm 0.008`$, is very close to the photometric value determined most recently by Lister et al. (2000), $`q_{ph}=0.53\pm 0.02`$, hopefully ending the long dispute about the $`q_{sp}`$ versus $`q_{ph}`$ disparity (Maceroni, 1986).

The ephemeris used to phase the observations was based on two sources: The period, as given in Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto., is from Nelson (2001) while the initial epoch, $`T_0=2,451,822.4267`$, was taken from Pribulla et al. (2001).

V523 Cas does not have a trigonometric parallax and was apparently too faint to be included in the Tycho–2 project. Thus, we refrain from a full discussion of its parameters, noting only that our spectral type, K4V, confirms previous classifications and agrees with $`(BV)=1.07`$ measured by Bradstreet (1981); in that work $`V_{max}=10.59`$. The (extrapolated) calibration of Rucinski & Duerbeck (1997, RD97) gives $`M_V=6.15`$ so that V523 Cas is one of the faintest known contact binaries. This gives a rather small distance to the binary, about 77 pc.

### 2.3 QW Gem

This contact system has been discovered by Hipparcos. It is somewhat faint for our telescope, so that the individual broadening functions are noisy, but the orbital solution is well defined (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.). This contact system is of the W-type when the Hipparcos ephemeris is used; this ephemeris gives a very small $`OC`$ deviations. The light curve has a moderately large amplitude of 0.41 mag.

The maximum light brightness, $`V_{max}=10.3\pm 0.1`$, is estimated from the transformed $`H_P`$ data and the Tycho–2 average data; this estimate is very approximate. The Hipparcos parallax has a large error, $`4.09\pm 5.42`$, which may indicate some undiscovered companions or just reflects the relative faintness of the star. Our spectral type, F8V, is consistent with the Tycho–2 average $`(BV)=0.48`$. $`M_V=3.36`$ determined from the very uncertain parallax is in agreement with $`M_V(cal)=3.55`$ derived using the RD97 calibration.

### 2.4 V921 Her

V921 Her is another Hipparcos discovery. This contact system is of the A-type with the more massive component eclipsed during the deeper minimum. The Hipparcos light curve is well defined, with an amplitude of 0.36 mag. The system is interesting in that it has a relatively long period of 0.877 days which is rather infrequent among contact binaries. The color index estimated from the Tycho–2 data, $`(BV)=0.27`$, points at a spectral type slightly later than our new classification, A7IV (the SIMBAD type is an even earlier A5), but there may be some interstellar reddening. The Hipparcos parallax, $`2.39\pm 0.97`$ mas and $`V_{max}=9.37`$, with an assumption of no reddening, give $`M_V=1.3\pm 0.9`$, which is consistent with the extrapolated RD97 calibration, $`M_V(cal)=1.1`$. This is one of the most luminous contact systems with $`P<1`$ day.

The components are only barely resolved in the broadening function which may indicate a low orbital inclination (a small $`\mathrm{sin}i`$). The orbit is moderately well defined (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.), partly due to the early spectral type and dearth of lines; both factors tend to increase the velocity errors. An interesting feature is a relatively large (in the absolute sense) mean systemic velocity of $`V_0=79`$ km s-1; the Hipparcos tangential velocities are small so that the large radial velocity is the main contributor to the spatial velocity of 81 km s-1.

### 2.5 V2357 Oph

This variable was discovered by Hipparcos. It was classified as a pulsating star with a period of 0.208 day which is two times shorter than the actual orbital period. The Hipparcos ephemeris was for the maximum light. Assuming that the Hipparcos initial epoch corresponded to the minimum just before the maximum used in the ephemeris ($`T_0=2,448,500.0451`$), the resulting $`OC`$ is quite small, see Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto..

The orbit is only moderately well defined (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.), partly because the system is a bit faint for our instrumentation, and partly because of the small separation of the components in the broadening functions. The small radial velocity semi-amplitudes $`K_i`$ and the small photometric amplitude of 0.12 mag. are all consistent with the low orbital inclination.

The Hipparcos parallax, $`5.37\pm 2.01`$ mas and $`V_{max}=10.43`$ lead to $`M_V=4.08\pm 0.82`$, while $`M_V(cal)=4.23`$ (RD97) for the average $`(BV)=0.80`$, as indicated by the Tycho–2 data. This color index is approximately consistent with our spectral type, G5V.

### 2.6 V1130 Tau

This bright ($`V_{max}=6.56`$) system was discovered by Hipparcos. It is the only non-contact system in this group of ten binaries. The broadening functions are exceptionally well defined and show (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.) two large, but detached components of similar size and brightness. The mass ratio resulting from our orbit (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.) is $`q=0.919\pm 0.008`$. The photometric variations observed by Hipparcos define a very nice light curve with rounded maxima, and an amplitude of of 0.38 mag. The distinction of which component is the hotter one and thus which of the minima should be called the primary minimum is not well established from the Hipparcos light curve. The broadening functions show the fainter component eclipsed during the Hipparcos’ primary minimum, which may indicate that the identification of the eclipses is actually incorrect.

The system begs to be analyzed in detail as it is one of the shortest period detached systems with the spectral type F3V. This spectral type is in agreement with the average $`(BV)=0.34`$, as derived from the Tycho–2 data. The Hipparcos parallax is relatively large (in fact the largest in this group of binaries) and well determined, $`15.35\pm 0.78`$ mas, leading to $`M_V=2.49\pm 0.12`$. The contact-binary calibration of RD97 obviously does not apply in this case; we can only note that the combined brightness corresponds to some 0.6 mag. above the Main Sequence, which simply reflects the detached binary nature of two similar stars.

We note that $`(M_1+M_2)\mathrm{sin}^3i=2.41\pm 0.03M_{}`$ for V1130 Tau is relatively large so that the orbital inclination must be close to 90 degrees.

### 2.7 HN UMa

Another Hipparcos discovery, HN UMa, is a rather commonly-looking contact system seen at moderately low orbital inclination because both the photometric variations (0.12 mag.) as well as the radial velocity semi-amplitudes, $`K_i`$, are rather small (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.).

The spectral type of F8V is somewhat early for the average $`(BV)=0.46`$ derived from the Tycho–2 data. The Hipparcos parallax, $`5.81\pm 1.39`$ mas, with the adopted $`V_{max}=9.80`$, give $`M_V=3.62\pm 0.53`$, while $`M_V(cal)=3.36`$ (RD97). The proper motion of the system (Tycho–2 data) is comparatively large, $`\mu _\alpha \mathrm{cos}(\delta )=+54.0\pm 1.9`$ mas and $`\mu _\delta =36.9\pm 1.8`$ mas, which together with $`V_0=37.1`$ km s-1 results in a relatively large spatial velocity of 65.0 km s-1.

### 2.8 HX UMa

The variability of HX UMa was discovered by Hipparcos. The light curve is of a typical W UMa-type system with an amplitude of 0.17 mag. Apparently, the ephemeris of Hipparcos predicts the photometric minima quite well with a very small $`OC`$ shift; the deeper one corresponds to the eclipse of the more massive star.

A close companion at the distance of 0.63 arcsec was discovered by Hipparcos. The re-determination by Fabricius & Makarov (2000) gave an average magnitude difference between the components of 3.00 mag. in the $`H_P`$ system; our determination of a weak third-light peak in the broadening function gives $`L_3/(L_1+L_2)=0.049\pm 0.004`$ or $`\mathrm{\Delta }m=3.27\pm 0.10`$ at the bandpass of our spectroscopy at 5184 Å. The signature of the third body was just barely detectable (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.) and, if not for the persistence through the orbital phases of the binary, it could be easily taken for noise in the broadening function. The radial velocity data for the third component are given in Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.. The average velocity, $`V_3=23.59\pm 0.77`$ km s-1, is slightly different from the center of mass velocity of the binary, $`V_0=19.88\pm 0.39`$ km s-1, but the difference may be of no significance as we cannot exclude any systematic bias in the measurements of $`V_3`$ on top of the much stronger feature of the close binary system. We do not see any “cross-talk” in the velocity determinations in the sense that $`V_3`$ does not depend on the orbital phase of the close binary (see Paper VII for a discussion of cases when such a dependence is observed). Our observations were made in two groups, centered on HJD 2,451,697.4 and 2,451,986.9 (see Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.). We do not see any significant change in $`V_3`$ between these epochs, with the average values of $`22.95\pm 1.12`$ and $`24.56\pm 1.00`$ km s-1.

The spectral type, F8V, is consistent with $`(BV)=0.44`$ derived from the Tycho–2 data. The Hipparcos parallax is poorly determined, $`6.68\pm 3.01`$, possibly reflecting the visual-binary nature of the system, so it cannot be meaningfully compared with the RD97 calibration which predicts $`M_V(cal)=3.32`$.

It was noted by Bartkevicius & Lazauskaite (1997) that HX UMa is not physically related to the nearby variable star DP UMa (HD 104513).

### 2.9 HD 93917

The W UMa-type variability of this star has been recently discovered by Lasala–Garcia (2001); we used the timing data from that work for the values of $`T_0`$ and $`P`$. The amplitude of light variations is 0.34 mag. so it is rather curious that the variability of this bright ($`V_{max}=9.02`$) star went undetected for such a long time. There are more curious circumstances around this star: It is one of very few HD stars not observed by the main Hipparcos project (but it was observed by Tycho–2). Also, there existed a large discrepancy between the $`(BV)=0.55`$ color index, as given in the SIMBAD database and confirmed by Tycho–2, and the SIMBAD spectral type of K0. We found that the spectral type is F9.5V which agrees very well with the color index. The predicted $`M_V(cal)=3.35`$ (RD97) cannot be checked because the Tycho parallax is poor, $`14.9\pm 11.4`$ mas.

Our observations define a very good spectroscopic orbit without any peculiarities (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.). The system belongs to the W-type of contact binaries, i.e. has a hotter less-massive component. There is absolutely nothing peculiar in the system except for the very fact that, being so bright, its variability remained undetected for such a long time. We note that Rucinski (2002c) (Sec. 8 & 12) estimated that even among brights stars within the interval $`7.5<V<8.5`$, perhaps as many as 45 contact systems still remain to be discovered.

### 2.10 NSV 223

The variability of NSV 223 was known since the discovery by Strohmeier et al. (1956). Recently, Verrot & Van Cauteren (2000) obtained a light curve showing a total, long-duration eclipse in the secondary minimum, which suggests that the mass ratio is small. Our radial velocity data (Figure Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.) fully confirm these predictions: The contact system is indeed of the A-type and the mass ratio is small, $`q=0.136\pm 0.011`$. The large semi-amplitudes $`K_i`$ suggest that the orbital inclination is close to 90 degrees. In spite of those properties promising a very good combined solution, the system is poorly known. What is known is $`V_{max}=10.86`$ from Verrot & Van Cauteren (2000); also our spectral type is F7V. The ephemeris given by the discoverers holds quite well (see Table Radial Velocity Studies of Close Binary Stars. VIII111Based on the data obtained at the David Dunlap Observatory, University of Toronto.).

## 3 SUMMARY

This paper presents radial-velocity data and orbital solutions for the seventh group of ten close binary systems observed at the David Dunlap Observatory. In this series, only V523 Cas has a history of previous radial-velocity studies. Seven system were discovered photometrically by the Hipparcos satellite mission while two binary systems, HD 93917 and NSV 223 have been identified as contact systems only very recently.

All systems of this group, with the exception of the detached binary V1130 Tau, are contact binaries of the W UMa type. All systems are double-lined (SB2) binaries with visible spectral lines of both components. V1130 Tau is interesting one because it is a very close yet detached consisting of two almost identical F3V stars on a tight orbit with $`P=0.799`$ days. In contrast, the early-type (A7) contact binary V921 Her has a relatively long period of $`P=0.877`$ days. A bright (9 mag.), recently discovered contact system HD 93917 is rather inconspicuous in its properties except for the fact that it remained undiscovered for such a long time; this confirms the prediction that many close binary systems still remain to be discovered among bright stars. V410 Aur and HX UMa have fainter companions detected spectroscopically.

Support from the Natural Sciences and Engineering Council of Canada to SMR and SWM is acknowledged with gratitude. The research has made use of the SIMBAD database, operated at the CDS, Strasbourg, France and accessible through the Canadian Astronomy Data Centre, which is operated by the Herzberg Institute of Astrophysics, National Research Council of Canada.

Captions to figures:


text:
# Computation of Lyapunov Exponents in General Relativity

## Abstract

Lyapunov exponents (LEs) are key indicators of chaos in dynamical systems. In general relativity the classical definition of LE meets difficulty because it is not coordinate invariant and spacetime coordinates lose their physical meaning as in Newtonian dynamics. We propose a new definition of relativistic LE and give its algorithm in any coordinate system, which represents the observed changing law of the space separation between two neighboring particles (an “observer” and a “neighbor”), and is truly coordinate invariant in a curved spacetime.

Chaos is a popular phenomenon in dynamical systems. One of its main features is the exponential sensitivity on small variations of initial conditions. The exhibition of chaos in the motion of Pluto makes it particularly attractive for scientists to investigate the dynamical behavior of the solar system.

In Newtonian mechanics, Lyapunov exponents (LEs), as a key index for measuring chaos in a dynamical system, can be calculated numerically with either the variational method or the two-particle method, and sometimes a mixture of the two. The former is more rigorous but one has to derive the variational equations of the dynamical equations of the system and integrate them numerically with the dynamical equations together. The latter is less cumbersome especially when one wants to compute the maximum LE only, which is the main index for chaos. We will discuss the extension of the two particle method to relativistic models in this letter.

Let $`𝐪(t)`$ and $`\dot{𝐪}(t)`$ be the position and velocity vector of a dynamical system. As a set of initial conditions is selected randomly and the corresponding trajectories are restricted in a compact region in the phase space, the classical definition of the maximum LE is

$$\lambda _N=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\frac{d(t)}{d(0)},$$
(1)
where the separation $`d`$ between two neighboring trajectories is required to be sufficiently small so that the deviation vector $`(\mathrm{\Delta }𝐪,\mathrm{\Delta }\dot{𝐪})`$ can be regarded as a good approximation of a tangent vector, and the distance $`d(t)`$ at time $`t`$ is of the form

$$d(t)=\sqrt{\mathrm{\Delta }𝐪(t)\mathrm{\Delta }𝐪(t)+\mathrm{\Delta }\dot{𝐪}(t)\mathrm{\Delta }\dot{𝐪}(t)}.$$
(2)
It is the Euclidian distance in the phase space. $`𝐪`$ and $`\dot{𝐪}`$ are in different dimension, so one has to carefully choose their units to assure that both terms in the expression of $`d(t)`$ are important. We suggest computing the LE in the configuration space rather than in the phase space, that is,

$`d^{}(t)`$ $`=`$ $`\sqrt{\mathrm{\Delta }𝐪(t)\mathrm{\Delta }𝐪(t)}`$ (3)
$`\lambda _N^{}`$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}\mathrm{ln}{\displaystyle \frac{d^{}(t)}{d^{}(0)}},`$ (4)
We argue that both $`\lambda _N`$ and $`\lambda _N^{}`$ are effective in detecting chaotic behavior of orbits since $`\mathrm{\Delta }\dot{𝐪}(t)=d/dt(\mathrm{\Delta }𝐪(t))`$ and they should have the same Lyapunov exponents.

In general relativity the above definition and algorithm are questionable. First, there is no unified time for all the reference systems. Secondly, the separation of time and space of the 4-dimensional spacetime is different for different observers. Furthermore, time and space coordinates usually play book-keeping only for events and are not necessary with a physical meaning. Consequently, one would get different values of $`\lambda _N`$ and $`\lambda _N^{}`$ in different coordinate systems.

Up to now the references for exploring chaos in relativistic models have been mainly interested in studying the dynamical behavior of black hole systems\[5-8\] and the mixmaster cosmology\[9-13\]. Most of them adopt the classical definition of LE. Here we will abandon the variational method in the case of general relativity for it is rather cumbersome to derive variational equations. Actually there exist general expressions for geodesic deviation equations but one has to derive the complicated curvature tensor. Furthermore, in many cases a particle does not follow a geodesic. For instance, a spinning particle in Schwarzschild spacetime doesn’t move along a geodesic. Therefore, we will concentrate on constructing a revised edition of the two-particle method. The classical algorithm of LE lacks invariance in general relativity, which depends on a coordinate gauge and even could bring spurious chaos in some coordinate systems. The chaotic behavior in the mixmaster cosmology\[9-13\] has been debated for decade or so. Using different time parameterization , one would get distinct values of $`\lambda _N`$. Especially, in a logarithmic time, chaos becomes hidden. Chaos, as an intrinsic nature of a given system, should not be affected by the choice of a coordinate gauge, and a chaos indicator, LE, should be defined as invariant under spacetime transformations. It means that the LE in general relativity should be expressed as a physical or so called proper quantity but not a coordinate quantity\[16-19\].

Now we consider a particle, called “observer”, moving along an orbit (not necessary to be a geodesic) in the spacetime with a metric $`ds^2=g_{\alpha \beta }dx^\alpha dx^\beta `$. In this letter Greek letters run from 0 to 3 and Latin letters from 1 to 3. The observer can determine if his motion would be chaotic by observing whether the proper distances from his neighbors are increasing exponentially or not with his proper time. Both the distances and time are observables and should not depend on the choice of a coordinate system. Then we can apply the theory of observation in general relativity to explore the dynamical behavior of the observer.

Fig. 1 illustrates the trajectories of the observer and its neighboring particle called “neighbor”. The initial condition of the neighbor is arbitrarily chosen as long as its 4-dimensional distance from the observer is small enough to assure that the neighbor is approximately in the tangent space of the observer. In an arbitrary spacetime coordinate system $`x^\alpha `$, one can derive the equations of motion of the observer and its neighbor. Here we have to notice that the independent variable of the equations has to be the coordinate time $`t`$ rather than their proper times $`\tau `$ because the two particles have different proper times but one unique independent variable has to be adopted when integrating numerically their equations of motion together. At the coordinate time $`t`$ the observer arrives at the point $`O`$ with the coordinate $`x^\alpha `$ and 4-velocity $`U^\alpha `$, and the corresponding proper time of the observer is $`\tau `$. At the same coordinate time $`t`$, its neighbor reaches the point $`P`$ with the coordinate $`y^\alpha `$ along another orbit. A displacement vector $`\delta x^\alpha =y^\alpha x^\alpha `$ from $`O`$ to $`P`$ should be projected into the local space of the observer. The space projection operator of the observer is constructed as $`h^{\alpha \beta }=g^{\alpha \beta }+c^2U^\alpha U^\beta `$ ($`c`$ represents the velocity of light). $`\stackrel{}{OP^{}}`$ represents the projected vector $`\delta x_{}^\alpha =h_\beta ^\alpha \delta x^\beta `$, and its length $`\stackrel{}{OP^{}}`$ is

$$\mathrm{\Delta }L(\tau )=\sqrt{g_{\alpha \beta }\delta x_{}^\alpha \delta x_{}^\beta }=\sqrt{h_{\alpha \beta }\delta x^\alpha \delta x^\beta }.$$
(5)
Here $`g_{\alpha \beta }`$ is calculated at $`x^\alpha `$. $`\mathrm{\Delta }L(\tau )`$ is the proper distance to the neighbor observed by the observer at his time $`\tau `$ and it is a scalar. Hence the maximum LE in general relativity is defined as

$$\lambda _R=\underset{\tau \mathrm{}}{lim}\frac{1}{\tau }\mathrm{ln}\frac{\mathrm{\Delta }L(\tau )}{\mathrm{\Delta }L(0)}.$$
(6)
Let $`\mathrm{\Sigma }_\tau `$ be the 3-dimensional subspace of the tangent space of the observer at $`O`$, which is orthogonal to the 4-velocity $`U^\alpha `$. It is the point $`P^{}`$ but not $`P`$ in $`\mathrm{\Sigma }_\tau `$ as long as $`P^{}`$ is close enough to the observer $`O`$. This tells us that $`\lambda _R`$ is coordinate invariant. The next is an argument for this point. Let us carry out a time transformation $`t\eta `$. For convenience, we express the projection operation as $`\stackrel{}{OP^{}}=H\stackrel{}{OP}`$, where $`H`$ is the space projection operator and $`H\stackrel{}{U}_O=0`$, where a subscript is added for $`\stackrel{}{U}`$ to describe the 4-velocity. Assume that the observer is located at the point $`O`$ at $`\eta `$ that corresponds to the old coordinate time $`t`$, then the neighbor at $`\eta `$ would situate at the point $`Q`$ but not necessary at $`P`$. Certainly we have to keep both $`P`$ and $`Q`$ inside an $`ϵ`$ neighborhood of $`O`$ (see Fig.1) and $`ϵ`$ is considered as a very small quantity. In this case we have $`\stackrel{}{U}_P=\stackrel{}{U}_O+O(ϵ)`$. Furthermore, we have $`\stackrel{}{OQ}=\stackrel{}{OP}\stackrel{}{QP}`$ and $`\stackrel{}{PQ}=\mathrm{\Delta }\tau _P\stackrel{}{U}_P+O(\mathrm{\Delta }\tau _P^2)=kϵ\stackrel{}{U}_P+O(ϵ^2)`$, where $`k`$ is a constant and the proper time interval of the neighbor between $`P`$ and $`Q`$, $`\mathrm{\Delta }\tau _P`$, is in the magnitude of $`ϵ`$. Hence, we get the relation $`H\stackrel{}{OQ}=H\stackrel{}{OP}+O(ϵ^2)`$. This shows that $`\lambda _R`$ is invariant with time transformations.

As far as the numerical implementation of the computation of $`\lambda _R`$ is concerned, the following notes are worth noticing. (1) The observer and the neighbor have different proper times, so a coordinate time $`t`$ should be adopted as the independent variable. The equations of the motion of the particles should be transformed to use $`t`$ as the independent variable. The variables to be computed step by step are the space coordinates and velocities (with respect to $`t`$) of the observer and its neighbor, and the proper time $`\tau `$ and $`d\tau /dt`$ of the observer. In total there are 14 variables. (2) As is well known, renormalization after a certain time interval is essential in this procedure to keep the distance between the observer and its neighbor small enough. The renormalization must be proceeded in the phase space though our LE is calculated in the configuration space. (3) One important difference during renormalization between the relativistic and classical cases must be noticed. Let $`\mathrm{\Sigma }_t`$ be the local 3-space in which all the points have the same coordinate time $`t`$. It is evident that the renormalization should be done in $`\mathrm{\Sigma }_t`$ otherwise the computation will commit an error. In Fig.1 $`P`$ is pulled back to $`M`$ and $`\stackrel{}{OM}=(\mathrm{\Delta }L(0)/\mathrm{\Delta }L(\tau ))\stackrel{}{OP}`$. The next integration step for the neighbor will start from the point $`M`$. It is obvious that $`\stackrel{}{OP}`$ is located in $`\mathrm{\Sigma }_t`$ because both points $`O`$ and $`P`$ are in $`\mathrm{\Sigma }_t`$. On the other hand $`\stackrel{}{OP^{}}`$ is in the 3-space $`\mathrm{\Sigma }_\tau `$ and a pull back from $`P^{}`$ to $`S`$ as renormalization is not correct. Our practice has proven this argument.

In order to check the validity of our scheme for calculating the relativistic LE, we are going to reexamine the core-shell system studied by Vieira & Letelier. In Schwarzschild coordinates $`(t,r,\theta ,\varphi )`$, the 4-metric for this system is of the form

$`ds^2`$ $`=`$ $`(1{\displaystyle \frac{2}{r}})e^Adt^2+e^{BA}[(1{\displaystyle \frac{2}{r}})^1dr^2+r^2d\theta ^2]`$ (7)
$`+e^Ar^2\mathrm{sin}^2\theta d\varphi ^2,`$
where $`A`$ and $`B`$ are functions of $`r`$ and $`\theta `$ only. Here we adopt nondimensional variables and take $`c=1`$. The metric does not explicitly depend on $`t`$ and $`\varphi `$ and has an energy constant $`E`$ and an angular momentum constant $`L`$, so test particles in free fall are actually in a system with two degrees of freedom with an integral $`U^\alpha U_\alpha =1`$. When $`A=B0`$, Eq.(7) represents the Schwarzschild spacetime in which test particles in free fall move in regular orbits due to integrability of the system. For simplicity, we discuss the dynamical behavior of geodesics in a black hole plus a dipolar shell, and let

$`A`$ $`=`$ $`2\sigma \mu \upsilon ,`$
$`B`$ $`=`$ $`\gamma _0+4\sigma \upsilon \sigma ^2[\mu ^2(1\upsilon ^2)+\upsilon ^2],`$ (8)
where $`\mu =r1`$ and $`\upsilon =\mathrm{cos}\theta `$. The Newtonian limit of the model resembles the Stark problem, which is fully integrable. However, as far as the relativistic model is concerned, Vieira & Letelier and Saa & Venegeroles have demonstrated strong chaos in Weyl coordinates using the Poincar$`\stackrel{´}{e}`$ surface of section, and Gu$`\stackrel{´}{e}`$ron & Letelier estimated its maximum LE $`\lambda _N=(3.2\pm 0.4)\times 10^4`$ with the classical definition of LE (see Eq.(1)) in prolate spheroidal coordinates.

We choose parameters as $`E=0.975`$, $`L=3.8`$, $`\sigma =2.5\times 10^4`$, and $`\gamma _0=\sigma ^2`$, and the initial conditions of the observer as $`r=32`$, $`\theta =\frac{\pi }{2}`$, $`\varphi =0`$, $`\dot{r}=0`$, and $`\dot{\theta }`$ from $`U^\alpha U_\alpha =1`$. As to its neighbor, an initial separation $`\mathrm{\Delta }r=10^8`$ is adopted, regarded as the best choice, and the others remain the same as the observer’s except $`\dot{\theta }`$. The values of $`E`$ and $`L`$ are carefully chosen to assure that the trajectories of the observer and its neighbors are bounded in a compact region. We integrate the geodesic equations using two integrators for comparison, Runge-Kutta -Fehlberg 7(8) and the 12th-order Adams-Cowell method, with a coordinate time step 0.01. When $`t`$ reaches $`10^6`$, we find $`\lambda _N=\lambda _N^{}=(2.2\pm 0.2)\times 10^4`$, which is close to the result of  as its integration time amounts to $`10^5`$. Meanwhile, we get $`\lambda _R=(2.8\pm 0.2)\times 10^4`$ by Eq.(6). One may notice that $`\lambda _N`$ and $`\lambda _R`$ are in the same magnitude. This is because the Euclidian distance $`d^{}(t)`$ and the invariant proper distance $`\mathrm{\Delta }L`$ differ not very much and so do the coordinate time $`t`$ and the proper time $`\tau `$. In fact, $`\tau `$ runs about $`9.174\times 10^5`$ when $`t`$ passes through $`10^6`$. A stronger gravitation by decreasing the angular momentum $`L`$ will increase the difference between $`\lambda _N`$ and $`\lambda _R`$, but the smallest stable circular orbit in Schwarzschild spacetime is $`r=6`$.

To verify the invariance of $`\lambda _R`$, we do a time transformation

$$t\eta =10t+r^2/2.$$
(9)
We obtain $`\lambda _N=(2.6\pm 0.2)\times 10^5`$ in the coordinate time $`\eta `$, while $`\lambda _R`$ retains the original value (see Fig.2 and Fig.3).

As a further experiment, we slightly change the model to put $`B/2A=2\sigma \mu \upsilon `$. This system is still chaotic in Schwarzschild coordinates $`(t,r,\theta ,\varphi )`$. To remove the singularity of this metric at the horizon, we go to Lemaître coordinates $`(T,R,\theta ,\varphi )`$ as

$`T`$ $`=`$ $`\kappa (t+2\sqrt{2r}+2\mathrm{ln}|{\displaystyle \frac{\sqrt{r}\sqrt{2}}{\sqrt{r}+\sqrt{2}}}|),`$
$`R`$ $`=`$ $`{\displaystyle \frac{2\kappa r}{3}}\sqrt{{\displaystyle \frac{r}{2}}}+T,`$ (10)
where the positive constant $`\kappa `$ may be chosen arbitrarily. Our numerical tests display that the larger $`\kappa `$ is, the smaller $`\lambda _N`$ becomes. Particularly, the Lemaître coordinates hide chaos for sufficient large $`\kappa `$ if $`\lambda _N`$ is adopted as a chaotic index. However, our $`\lambda _R`$ does not vary with the parameter $`\kappa `$.

We would like to emphasize again that the computation of $`\lambda _R`$ can be applied, whether the particles move along geodesics or not. The theory of observation in general relativity is not relevant to the 4-acceleration of the observer. For example, this method can be used to identify chaos in compact binary systems.

In this letter we concentrate on calculating an invariant maximum LE in general relativity, but this discussion can be easily extended to find all the relativistic LEs numerically with the technique proposed by Benettin et al. (also see ). They suggest choosing an initial set of orthonormal tangent vectors as a base of the tangent space of the observer, then compute the evolution of the volume determined by these vectors. To avoid two vectors getting close to each other under evolution they use Gram-Schmidt procedure to orthonormalize the base after each time step. The only change in a relativistic model is in computing the evolution of a tangent vector, which can be realized by the technique in this letter. As to the inner product in the orthonormalizing procedure and the norm computation, one has to use the Riemanian inner product in place of the Euclidean one.

We are very grateful to Dr. Xin-lian Luo of Nanjing University for his helpful discussion. This research is supported by the Natural Science Foundation of China under contract Nos. 10233020 and 10173007.