text:
# Disks surrounding Cataclysmic Binaries

## 1. Introduction

The standard picture for the evolution of Cataclysmic Variables (CV), now nearly 2 decades old, has been reasonably successful as a framework for interpreting the phenomenology of CV. In this interpretation, angular momentum loss by gravitational waves and magnetic braking causes cataclysmic binaries to evolve from longer to shorter orbital periods ($`P`$), with magnetic braking dominating above the ‘period gap’ ($`P>3`$ hr), and a lower braking rate below this gap ($`P<2`$ hr). This picture explains some overall statistical properties, for example the existence of the period gap itself and the higher mass transfer rate above the gap compared to short period systems.

An additional element of the phenomenology, independent of the magnetic braking scenario but not in conflict with it either, is the existence of outflows from the accretion disk. Evidence for outflows comes from the profiles in UV-lines seen in outbursts of Dwarf Novae (e.g. Woods et al. 1992), possibly the single-peaked line profiles of the novalike (NL) and SW Sex binaries (see however Hellier 2000 for an alternative explantation of these line profiles), and other indications for outflow (e.g. Long et al. 1994).

## 2. Anomalies

A number of observations, however, are hard to fit into the standard picture:

* 1 A major problem is the very large spread in inferred mas transfer rates at a given orbital period (Patterson 1984; Warner 1987). The braking scenario predicts a close relation between mass transfer rate and orbital period.
* 2 An extreme form of the spread in mass transfer rate is shown by the short-period super soft sources (SSS). Though some SSS have long orbital periods and their large transfer rates can be explained as the result of expansion of the secondary by nuclear evolution, several have short orbital periods, in the same range as the novalike systems. The secondaries at such periods are too small to have significantly evolved. An unknown process causes these binaries to transfer mass at rates sufficient ($`\mathrm{3\hspace{0.17em}10}^7\mathrm{M}_{}\mathrm{yr}^1`$) to sustain the observed steady burning of Hydrogen at the surface of the White Dwarf primary star.
* 3 A second unexpected fact is the distribution of the novalike systems, which cluster at periods (around 4-6 hr) just above the period gap, whereas the braking scenario predicts them to predominate at longer periods, where the predicted mass transfer rate are highest.
* 4 The single peaked profiles of the emission lines in the SSS, the SW Sex stars, the AM CVn stars and many NL systems show that in some cases the emission lines are not predominantly formed in the accretion disk, but in optically thin material of large emission measure somewhere else in the system. The phenomenon is particularly striking in some AM CVn systems, (cf. Warner 1995), where a narrow emission line is seen on top of a broad double peaked emission line from the accretion disk. This phenomenon has been interpreted in terms of circumbinary material already by Solheim and Sion (1994). A similar phenomenon may be the low-velocity emission peaks in supersoft sources like QR And (Deufel et al. 1999).
* 5 ‘Additional light’ in systems with parameters such that mid-eclipse of the primary star and accretion disk should be total. In some of these systems the light remaining at mid-eclipse can not be explained by the spectral type of the secondary star. In the case of UU Aqr the spectrum of the additional light has been determined; it shows the characteristics of optically thin emission (Baptista et al. 1998).
* 6 Comparison of the spectral types of the secondary stars with theoretical models (Baraffe and Kolb 2000) shows them to be significantly cooler, on average, than main sequence stars, as if they were less massive and overexpanded as a result of a recent prolonged period of large mass transfer.
* 7 Superhumps above the period gap. The so-called superhumps, a feature in the lightcurve that drifts in orbital phase, is thought to be due to a resonance in the accretion disk. This resonance can only occur at low mass ratios, such as occur in the short-period systems. They are also seen in several novalike systems, however (Patterson, 1999). This would be another indication that the mass of the secondary in such systems above the period gap is often smaller than expected from the standard braking scenario.
* 8 Pile up at the minimum period. The standard braking scenario predicts that the orbital evolution of CVs slows down as they approach the minimum orbital period, after which their orbital period increases again at a very slow rate. Most old CVs are predicted to end up somewhere near this minimum period. The observed period distribution does not show such a spike in the distribution (cf. Kolb & Baraffe 1999). Though the observations are probably biased against detection of these systems because of their low quiescent luminosity and the long recurrence time of their outbursts, the discrepancy is large, and can probably not be explained entirely by these selection effects (Kolb 2001; Marsh et al. 2001; and references therein; Gänsicke et al. 2001).
* 9 Value of the minimum period. The theoretically predicted value of the minimum orbital period of CVs (excepting thiose with helium star companions) is around P=70 min, the observed value around 80 min. Stellar models are believed to be sufficiently accurate, but agreement with observations would be obtained if the mean mass transfer rate of the short-period systems were at least a factor 3 higher than assumed in the standard braking scenario.
* 10 The rate of magnetic braking. The standard formulas used for angular momentum loss by magnetic stellar winds is somewhat of an extrapolation, based on loss rates inferred from stars rotating much more slowly than the secondaries of CVs. Recent work by Andronov, Pinsonneault and Sills (2001) indicates that this extrapolation overestimates the angular momentum loss rate in CVs by as much as 2 orders of magnitude. If this is correct, the observed mass transfer rates in CVs can not be due to magnetic braking.

## 3. Interpretation

Though some of the anomalies listed are possibly red herrings which may find unrelated explanations, they all fit qualitatively into the suggestion made here, a CB disk, i.e. a disk surrounding the binary. If these exist, they would be a major, thus far unrecognized component of Cataclysmic binary systems.

A CB disk would be an additional source of angular momentum loss for the binary. The tides raised in the CB disk by the orbiting binary transfers angular momentum to it, and prevent it from spreading inward. Depending on the surface density of the CB disk, the angular momentum loss from the binary would differ, and with it the mass transfer rate. The SSS would, in this interpretation, be the relatives (Greiner et al. 2000) of NL and SW Sex systems, the difference being the presence of a more massive CB disk.

The increased mass transfer rates caused by CB disks would be systematic, since a CB disk would dissipate only very slowly by spreading. This might explain the apparent overexpansion of the secondary stars. A higher mass transfer rate would also help to explain the discrepancy between the observed and predicted values of the period minimum.

The obvious question is how a CB disk would be formed. One possibility would be that it is a ‘fossil’: a remnant from a late stage of the common envelope (CE) evolution phase in which the cataclysmic binary was formed (for a review see Taam & Sandquist 2000). This could happen as a ‘fall back’ process at the time the common envelope becomes optically thin. Given the difficulty of computing the CE evolution, it is hard to estimate the likelihood of such a formation channel for CB disks.

A second possibility is that CB disks form as a by-product of mass loss from the system. The single peaked line profiles in SW Sex stars indicate outflow velocities up to 2000 km/s, which is well above the escape speed from the binary. There is large optical depth in these lines at lower velocities as well, however, so that a significant outflow component at velocities of the order of the escape speed from the binary (a few 100 km/s) or lower may also exist. This would be material that would ‘hang around’ the binary rather than being ejected. It would naturally form a CB disk by cooling and settling to the equatorial plane. It is hard to quantify from the observations how much material might condense into a CB disk in this way. Instead, we have used model calculations in which this uncertainty is parametrized. These show that a relatively minor amount would suffice to build up substantial CB disk.

## 4. Model calculations

For the model calculations (Spruit and Taam 2001, Taam and Spruit 2001, Dubus et al. 2001) we have made use of the observation that evidence for outflows is strongest in systems with high mass transfer rates: dwarf novae in outburst, the SW Sex stars, supersoft sources and the AM CVn stars. This suggests that mass loss by outflows increases with mass transfer rate. We have made the simple assumption that the two are proportional. In addition, we have assumed that of the outflow, a fixed fraction settles into a CB disk. The mass input rate $`\dot{M}_{\mathrm{CB}}`$ into the CB disk is thus a given fraction $`\delta `$ of the mass transfer rate $`\dot{M}`$.

### 4.1. Feedback

The consequence of the assumed proportionality between mass input into the CB disk and the mass transfer rate is a positive feedback: once a CB disk is present, the angular momentum it extracts from the binary enhances the mass transfer rate from the secondary to the primary. This increases the mass input rate into the CB disk, increasing its surface density. With this higher density the increasing torque on the binary in turn increases the mass transfer rate. Though the mass transfer might initially be caused by the standard magnetic braking process alone, the mass building up in the CB disk would eventually make this feedback strong enough to become unstable. From this point on, the mass transfer rate and the input rate into the CB disk would increase in step exponentially. In this way the large mass transfer rates in for example the short-period supersoft sources could develop. In our interpretation, these would be descendants of novalike systems, as suggested before by Greiner et al. (2000).

We call the phase of high mass transfer resulting from this instability ‘accelerated mass transfer’ induced by the CB disk. The calculations show that the orbital period typically ‘bounces’ shortly after accelerated mass transfer sets in (if the secondary is not significantly evolved) . The more rapid mass loss from the secondary causes it to expand, forcing the orbital period to increase. This contrasts with the evolution under standard braking, which results in an essentially monotonically decreasing orbital period.

Whether the effect of a CB disk on the binary becomes significant depends mostly on the mass input rate into it. Analytic results based on the above assumptions (Spruit and Taam 2001) suggested that a minimum value of $`\delta `$ around $`10^3`$ would be sufficient to cause accelerated mass transfer to occur. The precise value is unfortunately somewhat sensitive to the physics determining the viscous spreading of the disk. Remarkably, the uncertain value of the viscosity itself, as parametrized through a constant $`\alpha `$-viscosity, has only a weak influence. This is because the torque exerted by the disk is proportional to $`\alpha \mathrm{\Sigma }_\mathrm{i}`$, where $`\mathrm{\Sigma }_\mathrm{i}`$ is the surface density at the inner edge of the disk. A larger value of $`\alpha `$ is compensated largely by the more rapid spreading of the disk, which lowers $`\mathrm{\Sigma }_\mathrm{i}`$. More important than $`\alpha `$ itself is the relative variation of the viscosity through the disk. Calculations using more detailed physics for the CB disk (Taam and Spruit 2001) indicated that a rather high value of $`\delta `$, around $`10^2`$, would be needed for accelerated mass transfer, but these did not include the effect of convection on the structure of more massive CB disks. This was included by Dubus et al. (2001), together with a more sophisticated method for calculating the viscous evolution of the CB disk. In these so far most realistic calculations values of $`\delta `$ around $`10^4`$ were found to be sufficient to cause a rapid acceleration to high transfer rates to happen within a Hubble time.

## 5. Consequences of the presence of CB disks

### 5.1. Dissolution of the secondary

In the standard magnetic plus gravitational radiation braking scenario, the mass transfer rate declines rapidly with decreasing secondary mass, such that complete transfer of the secondary mass does not happen within a Hubble time. With torques due to a CB disk the evolution of the secondary is more interesting. In the ‘accelerated’ transfer phase, the angular momentum loss is determined by the surface density of the CB disk, rather than by the secondary star. Hence the mass transfer rate stays finite as the secondary mass vanishes. As a consequence, the secondary dissolves completely within a finite, relatively short time once a phase of accelerated mass transfer has set in, and leaves a single white dwarf.

This would have significant effects on the distribution of systems with orbital period. If a substantial fraction of all systems go through a CB-induced phase of accelerated mass transfer, they would disappear from the CV population, and the pile-up of systems towards the period minimum around 80 min. would be replaced by a gradual decline of the number of systems. This would be more in line with current observations (cf Kolb 2001).

### 5.2. White dwarfs with planets?

After dissolution of the companion, the white dwarf primary would be left as a ‘white widow’, with a CB disk still surrounding it. The mass in this disk could be substantial, perhaps as much as $`10^310^2\mathrm{M}_{}`$. As this disk cools, it may well form planets in the same way as a protoplanetary disk. It remains to be seen whether such planets would be detectable, for example through Doppler shifts of the NLTE cores of the Balmer lines of the WD.

### 5.3. The period gap

The CB interpretation given here implies a large systematic spread in mass transfer rates. \[As opposed to cyclic mass transfer variations such as proposed by King et al. (1995). In this mechanism the average mass transfer rate does not spread much\]. Though this spread explains a number of facts such as the observed varying degree of overexpansion of the secondaries, it is not easily compatible with the explanation of the period gap in terms of disrupted magnetic braking. The period gap gets ‘washed out’ if the long-term transfer rate differs too much between systems. It remains to be seen whether a plausible alternative explanation of the gap can be found within the CB scenario.

## 6. Observability

One may wonder how it is possible that something as serious as circumbinary disks around CVs could have escaped detection so far. On closer inspection, however, it turns out that such disks would in fact be rather difficult to observe. On the one hand, the luminosity of a CB disk is low. Though it extracts a lot of angular momentum, this happens at a large distance in a shallow gravitational potential, and is thus accompanied by little energy dissipation. This makes them inconspicuous against the very bright accretion disk, which processes material deep in the potential well of the primary. Secondly, the size of the CB disk is such that its spectral energy distribution peaks in the poorly accessible mid IR (at 3-10$`\mu `$m, cf. the models of Dubus et al., 2001). Finally, model calculations show that the CB disk probably lies in the shadow of the outer part of the accretion disk. Enhanced brightness through illumination of the disk by the bright central regions of the accretion disk is thus also unlikely.

## References

Andronov, N., Pinsonneault, M., & Sills, A. 2001, Submitted to ApJ (Astro-ph/0104265)

Baptista, R., Silveira, C., Catalan, S., Steiner, J.E., & Horne, K. 1998, IX Latin American Regional IAU Meeting, ”Focal Points in Latin American Astronomy”, Eds: Aguilar, A.; Carraminana, A.; Revista Mexicana de Astronomia y Astrofisica Serie de Conferencias., E26

Baraffe, I., & Kolb, U. 2000, MNRAS, 318, 354

Deufel, B., Barwig, H., Simic, D., Wolf, S., & Drory, N. 1999, A&A, 343, 455

Dubus, G., Taam, R.E. & Spruit, H.C. 2001, ApJ, submitted

Gänsicke B., Hagen, H.-J. & Engels, D. 2001, astro-ph/0100375

Greiner, J., Orio, M., & Schwarz, R. 2000, A&A, 355, 1041

Hellier, C. 2000, NewAR, 44, 131

Harrison, T.E., McNamara, B.J., Szkody, P. & Gilliland, R.L. 2000, AJ, 120, 2649

King, A.R., Frank, J., Kolb, U. & Ritter, H. 1995, ApJ, 444, L37

Kolb, U. 2001, astro-ph/0110253

Kolb, U. & Baraffe, I. 1999, MNRAS, 309, 1034

Long, K.S., Wade, R.A., Blair, W.P., Davidsen, A.F., & Hubeny, I. 1994, ApJ, 426, 704

Marsh, T.R., Morales-Rueda, L., Steeghs D., Maxted, P., Kolb U., Boyle B., Croom, S., Loaring N., Miller L., Outram, P. Shanks, & T. Smith, R. 2001, astro-ph/0108334

Patterson, J. 1984, ApJS, 54, 443

Patterson, J. 1999, in Disk Instabilities in Close Binary Systems, Disk Instabilities in Close Binary Systems., Eds. S. Mineshige and J.C. Wheeler, Frontiers Science Series No. 26. Universal Academy Press, p.61

Solheim, S.-E., & Sion, E. M. 1994, A&A, 287, 503

Spruit, H.C., & Taam, R.E. 2001, ApJ, 548, 900

Taam, R.E. & Sandquist, E.L. 2000, ARA&A, 38, 113

Taam, R.E., & Spruit, H.C. 2001, ApJ, 561, 329

Warner, B. 1987, MNRAS, 227, 23

Warner, B. 1995, Cataclysmic Variable Stars, CUP, Cambridge, Ch. 9.7.2

Woods, J.A., Verbunt, F., Collier Cameron, A., Drew, J. E., & Piters, A. 1992, MNRAS, 255, 237


format: text

text:
# Going with the Flow: a Lagrangian approach to self-similar dynamics and its consequences

## 1 Introduction

When one tries to computationally locate constant shape traveling wave solutions in an equation that supports them, the standard procedure is to go to a traveling frame, so as to change the time-dependent problem into a steady state one. Since, however, the correct speed is not, in general, known a priori, two possibilities exist: either we will see no steady states (having guessed the wrong speed), or we will see infinitely many (if we were fortunate in guessing the correct speed). The simple solution to writing fixed point algorithms for this problem is to realize that one is allowed to pick one out of the the one-parameter infinity of steady state solutions of the PDE in the traveling frame. One-parameter here corresponds to traveling (translational invariance, a continuous symmetry group) in one spatial dimension; the concept generalizes directly to traveling in more than one spatial dimensions. Computationally one augments the system by writing an additional (possibly nonlocal) condition, called a “pinning condition” . We then have enough equations to solve for the one solution and the correct traveling speed simultaneously.

Rowley and Marsden  introduced recently a template-based technique for “factoring out” translational invariance; their work was motivated by the so-called Karhunen-Loeve expansion in the case of PDEs with traveling wave solutions. This was done by choosing a (more or less arbitrary) function, the template, and using it to systematically pick out one out of the one-parameter infinity of elements of the group orbit in a manner that will be explained below. We have recently implemented a discrete time version of this template-based approach to perform what we call “effective” bifurcation analysis of evolution equations in complicated media .

An important contribution of the work of Rowley and Marsden is that through the template method they allow the pinning to be done in a dynamic way. This means that we do not just implement some computational contraction mapping, and then converge on its fixed point, which is the solution traveling with constant speed and constant shape. Through the template we actually have a way to “factor out” the translational invariance in continuous time, and thus get an equation in an appropriate frame at all times. This is quite interesting since, often, what we want to find is not just the final solution, but we also want a view –“unobstructed” by the symmetry group– of the physical dynamics of the approach to the traveling wave. It is of course, in principle, possible to first find the traveling solution by some independent means (e.g. analytically) and then factor it out from the dynamics. The template approach allows us to do the “factoring out” of the group orbit naturally and continuously in time, without a priori knowledge of either the traveling solution or its speed. Moreover, as a simple byproduct it yields at the end the solution and the correct speed.

In this paper we extend the template approach to another continuous symmetry group: self-similarity. We take advantage of the fact that using a template to factor out self-similarity, we actually obtain a dynamic evolution equation. This equation not only gives the self-similar solutions, but it actually gives a view of the dynamics of approach to the self-similar solution with the self-similarity conveniently factored out.

The remainder of the paper is organized as follows: in the following section we show how to derive the transformed evolution equation for the scaled variable. We then illustrate the procedure through the computation of two well-known solutions (the Barenblatt solution, a self-similar solution of the first kind, and the Graveleau solution, a self-similar solution of the second kind) through the dynamics of the transformed equation. We conclude by a discussion of what we perceive as the implications of the approach, including alternative (to integration) computational approaches to the location of self-similar solutions. We also discuss briefly the discrete-time implementation of the approach on the original (as opposed to the transformed) equation. This has possibly important implications because of the recently proposed “coarse” integration/bifurcation analysis techniques for analyzing microscopic (e.g. Molecular Dynamics (MD) or Monte Carlo (MC)) timesteppers. It is conceivable that the computer-assisted analysis of self-similarity proposed here can be carried out even in cases where coarse, macroscopic equations describing the evolution of moments of molecular distributions conceptually exist but are not available in closed form.

We also note that what we propose here shares many elements with the approaches described, for example, in page 326 of the book of Goldenfeld , and in Chapter 6 of , including references to the work of the authors of the book, as well as Papanicolaou, Zakharov and their groups. We are currently exploring the similarities and differences of our work with that literature (see also ).

Because of the long ties, over the years, of two of the authors with the University of Minnesota, and because this work came to fruition in an office at U of M., we will call the transformed dynamics developed here the “MN-dynamics”.

## 2 Self-similar solutions

Self-similar solutions play an important role in the development of the theory of nonlinear evolutions equations. In addition to providing exact and sometimes even explicit solutions they often describe the asymptotic behavior at large times, or the asymptotic form of solutions in the neighborhood of some important change in behavior such as focusing or blow-up at some finite critical time. Often a self-similar solution to an autonomous evolution equation in the variables $`(x_1,x_2,\mathrm{},x_n,t)`$ is a function of the form

$$\left|tt^{}\right|^\beta F(\frac{x_1}{\left|tt^{}\right|^{\alpha _1}},\mathrm{},\frac{x_n}{\left|tt^{}\right|^{\alpha _n}}),$$
(1)
where the similarity exponents $`\alpha _1,\mathrm{},\alpha _n,\beta `$ and the function $`F`$ must be determined from the equation together with appropriate boundary and initial conditions. Here $`t^{}`$ is either the known initial time in which case we consider $`t>t^{}`$, or the unknown critical time so that we consider $`t<t^{}`$. In view of the time translation invariance we can always take $`t^{}=0`$. In some problems the exponents can be obtained a priori from scaling arguments and conservation laws. This is referred to as self-similarity of the first kind. However, it is often the case that the exponents cannot be gotten a priori and are usually obtained by solving what amounts to a nonlinear eigenvalue problem for the function $`F`$. This is the case of self-similarity of the second kind. A cogent account of the theory of self-similar solutions with many illuminating examples can be found in Barenblatt’s book .

In the next section we will derive the MN-dynamics equations for the construction of self-similar solutions. In section 4 we illustrate the method by applying it to the computation of the Barenblatt solution (first kind self-similarity) and the Graveleau solution (second kind self-similarity) of the porous medium equation (pme)

$$\frac{u}{t}=\mathrm{\Delta }(u^m)$$

where $`m>1`$ is a constant. For properties of the porous medium equation see .

## 3 MN-Dynamics

Consider the partial differential equation

$$u_t=D_x\left(u\right),$$
(2)
where the generally nonlinear operator $`D`$ acts on the variable $`x`$. We restrict ourselves to operators for which there exist two constants $`a`$ and $`b`$ such that for all positive $`A`$ and $`B`$

$$D_x\left(Bf\left(\frac{x}{A}\right)\right)=A^aB^bD_y\left(f\left(y\right)\right)\text{, where }y=\frac{x}{A}.$$
(3)
For these operators the PDE may have a one-parameter family of self-similar solutions

$$u(x,t)=s^\beta U\left(\frac{x}{Cs^\alpha }\right),$$
(4)
where $`C>0`$ is the parameter, $`s=\left|tt^{}\right|`$, and $`U`$ satisfies the ODE

$$\sigma C^a\left(\beta U\alpha yU_y\right)=D_y\left(U\right)\text{ in the variable }y=\frac{x}{Cs^\alpha }$$
(5)
with $`\sigma =sgn(tt^{})`$. The exponents $`\alpha `$ and $`\beta `$ satisfy the scaling condition

$$\beta 1=a\alpha +b\beta .$$
(6)
We now present an alternative method of computing similarity solutions, which factors out the similarity exponent $`\alpha `$. We start with the general scaling

$$u(x,t)=B(s)w(\frac{x}{A(s)},\tau (s)),$$
(7)
where $`A,B,`$ and $`\tau `$ are unknown functions. Using Eq. (3) we obtain the PDE for $`w`$

$$\sigma \left(\frac{B_s}{B}w\frac{A_s}{A}yw_y+\tau _sw_\tau \right)=A^aB^{b1}D_y(w)$$
(8)
Given an appropriate template function $`T(y)`$ we want to adjust $`A(s)`$ so that $`w`$ satisfies the orthogonality condition

$$_{\mathrm{}}^+\mathrm{}w(y,\tau )T(y)𝑑y=0.$$
(9)
Multiplying Eq. (8) by $`T(y)`$ and integrating we obtain

$$\frac{A_s}{A}=\sigma A^aB^{b1}I(\tau ),$$
(10)
where

$$I(\tau )=\frac{_{\mathrm{}}^+\mathrm{}D_y(w(y,\tau ))T(y)𝑑y}{_{\mathrm{}}^+\mathrm{}yw_y(y,\tau )T(y)𝑑y}.$$
(11)
To determine $`B`$ we need an additional pinning condition. For numerical convenience we impose the local condition $`w(p,\tau )=1`$ at some fixed point $`x=p`$. For example, $`p`$ may be a symmetry point where $`w_y(p,\tau )=0`$ for all $`\tau `$. Evaluating Eq. (8) at $`y=p`$ and substituting from Eq. (10) we obtain

$$\frac{B_s}{B}=\sigma A^aB^{b1}\left(D_y(w(p,\tau ))pw_y(p,\tau )I(\tau )\right).$$
(12)
Substituting Eqs. (10) and (12) in Eq. (8) we obtain

$$w_\tau +I(\tau )\left(yw_ywpw_y(p,\tau )\right)+wD_y(w(p,\tau ))=D_y(w)$$
(13)
provided that

$$\tau _s(s)=\sigma A^a(s)B^{b1}(s).$$
(14)
In view of Eq. (14) we can rewrite $`A`$ and $`B`$ as functions of $`\tau `$ instead of $`s.`$ Let $`\stackrel{~}{A}(\tau )=A(s(\tau ))`$ and $`\stackrel{~}{B}(\tau )=B(s(\tau ))`$. Then it follows from Eqs. (10), (12), and (14) that

$$\frac{\stackrel{~}{A_\tau }}{\stackrel{~}{A}}=I(\tau )\text{ and }\frac{\stackrel{~}{B_\tau }}{\stackrel{~}{B}}=D_y(w(p,\tau ))pw_y(p,\tau )I(\tau ).$$
(15)
Moreover, the physical time $`s=\left|tt^{}\right|`$ is given as a function of $`\tau `$ by

$$\frac{ds}{d\tau }=\sigma \stackrel{~}{A}^a(\tau )\stackrel{~}{B}^{1b}(\tau ).$$
(16)
Starting from a suitable initial condition $`w(y,0)`$, we integrate Eq. (13) for $`\tau \mathrm{}`$. If $`w(y,\tau )`$ converges to a steady state $`W(y)`$, there is a similarity solution to Eq. (2). In particular, $`W`$ satisfies

$$k\left(yW_ypW_y(p)W\right)+WD_y(W(p))=D_y(W),$$
(17)
where

$$k=\underset{\tau \mathrm{}}{lim}I(\tau )$$
(18)
as well as Eq. (5). Comparing coefficients in Eqs. (5) and (17) we find that

$$\frac{\alpha }{\beta }=\frac{k}{D_yW(p)kpW_y(p)}.$$
(19)
Thus the similarity exponents $`\alpha `$ and $`\beta `$ are determined by the linear system consisting of Eqs. (6) and (19).

In problems where there is a finite critical time $`t^{}`$ we are interested in $`t<t^{}`$ so that $`\sigma =1`$ and $`s=t^{}t`$. We integrate Eqs. (15) and (16) to obtain

$$t^{}=t(\tau =0)+\underset{\tau \mathrm{}}{lim}_0^\tau \stackrel{~}{A}^a(\eta )\stackrel{~}{B}^{1b}(\eta )𝑑\eta .$$
(20)
A remarkable property of this method is that is able to capture both self-similar solutions of the first and second kind, as we show in the next section.

We note that the new PDE Eq. (13) can also be obtained by repeatedly evolving Eq. (2) an infinitesimal time interval $`dt`$, and then rescaling the $`x`$ and $`u`$ variables such that they satisfy $`u(L,t)=1`$ and Eq. (9).

## 4 Examples

Here we give two examples of numerical integration of Eq. (13) to obtain self-similar solutions to the porous medium equation with $`m=2`$.

The 1D Barenblatt solution is the self-similar solution of the first kind to the porous medium equation whose initial datum is a Dirac mass at the origin. This solution is known explicitly  and

$$\alpha =\beta =1/3.$$

The operator in this case is $`D_x(u)=(u^2)_{xx}`$ and the scaling condition Eq. (6) reduces to $`\beta =2\alpha 1`$. We discretized the equation with centered finite differences. The time integration is carried out with an explicit Euler scheme on a domain $`[0,3]`$ discretized with $`600`$ points. From the linear stability analysis of the numerical scheme, the variable time step must satisfy $`\mathrm{\Delta }t<0.5/|w(y,\tau )|`$ on the domain of integration. To compute the Barenblatt solution we choose the template function

$$T(y)=\{\begin{array}{cc}1& \text{for }|y|1\\ 1& \text{for }|y|>1\end{array}$$

The orthogonality condition Eq. (9) means that the dimensionless mass within $`|y|1`$ equals the mass outside this interval. For the initial datum we take

$$w(y,0)=\{\begin{array}{cc}1& \text{for }|y|2\\ 0& \text{for }|y|>2\end{array}$$

which satisfies the orthogonality condition. Fig. 1 shows the convergence of the numerical solution to the steady state. The numerical value of the similarity exponent given by Eq. (19) is $`\alpha =0.333371`$ for $`\tau =10`$.

The Graveleau solutions to the porous medium equation form a 1-parameter family of axi-symmetric focusing solutions which are self-similar of the second kind . To construct them by integrating the MN-dynamic equation we write the equation in $`(r,t)`$ coordinates. For axi-symmetric diffusion the operator is $`D_r(u)=(u^2)_{rr}+(u^2)_r/r`$ and the scaling condition is again $`\beta =2\alpha 1`$. We approximate the solution that has an infinite support with a numerical solution in a finite interval $`[0,10]`$, discretized with $`2000`$ gridpoints. We use the boundary conditions $`w(0,\tau )=0,w_y(L,\tau )=0`$, and in Eq. (13) we set $`L=10`$.

The template function is chosen to be

$$T(y)=\{\begin{array}{cc}1& \text{for }0|y|7\\ 1& \text{for }7<|y|10\end{array}$$

and the initial datum is

$$w(y,0)=\{\begin{array}{cc}0& \text{for }0|y|4\\ 1& \text{for }4<|y|10\end{array}$$

which is orthogonal to $`T(y)`$ on the integration domain. We show the evolution of $`w`$ in Fig. 2. As $`\tau \mathrm{}`$ the solution tends to a steady state. In this case the corresponding physical time tends to a finite value $`t^{}`$ which is the time at which the solution to the original initial value problem first becomes positive in the neighborhood of the origin. $`t^{}`$ is called the focusing time. To determine the similarity exponent we use Eqs. (6) and (19). For 2 space dimensions the correct value of $`\alpha =0.856333\mathrm{}`$ and the MN-dynamics gives the approximate $`\alpha =0.85695`$.

The accuracy would be improved by using more gridpoints and enlarging the domain.

## 5 Discussion

Locating the self-similar solutions and their exponents by integration of the MN-equation is, of course useful. It is, however, only one of the ways that we can now find these solutions. Now that the self-similar solutions have become steady states of the MN-equation, we can bring all the tools of nonlinear dynamics to bear on their study. In particular, we can find self-similar solutions not only through the fact that they may be attracting (by integration), but also through contraction mappings like Newton’s method (this would be interesting in physical problems like the Landau-Gizburg equation ). Indeed, we can do bifurcation / continuation / stability calculations of the closed non-local PDE and find all types of stable, unstable, bifurcating, possibly even more complicated (e.g. periodic in the transformed frame, or even more complicated, chaotic in the transformed frame) self-similar solutions. They are going to be elements of the global attractor of the MN-equation, or, in general, invariant objects for the MN-dynamics. For example, in the traveling case, a limit cycle in the transformed equation would correspond to a modulated traveling wave — so, possibly, a limit cycle of the transformed equation would correspond to a ”periodically self-similar” solution for the self-similar problem. These might arise in Hopf bifurcations of the MN-equation, and might be associated with complex exponents; we are not sure at the moment what all these additional invariant objects might be, but it will be interesting to study them. We are certainly aware of cases that should correspond to steady state bifurcations in the MN-dynamics.

More importantly, however, the Jacobian of the linearization of the MN-equation at its steady states or limit cycles, should it prove to be a well defined object, will tell us what the stability in the self-similar frame is, with the self-similarity factored out. It will help us understand asymptotic rates of approach, bifurcations of new self-similar solutions, possibly find basins of attraction of different self-similar solutions as stable manifolds of saddle-self-similar solutions. In particular, it is conceivable that in the case of Type-II self-similar solutions, continuation and the use of this lenearization might facilitate the study of the post-focussing regime. There are many interesting singularities whose “other side” we would like to explore. The porous medium equation is one of them; possible examples from cosmology  also come to mind.

We conclude with a few (speculative at the moment) comments.

First, when a successful transformation turns the problem of finding self-similar solutions into finding steady states, limit cycles, chaotic attractors, or in general, elements of the global attractor of the transformed dynamics, we can bring all the machinery of dynamical systems / bifurcation theory to bear on computing self-similar solutions and their transitions (bifurcations in the tranformed frame). In particular, techniques much more powerful than simple integration (e.g. contraction mappings like the Newton method, continuation, and multiparameter bifurcation theory) can be brought to bear on the transformed problem.

Second, in everything we have done so far, we first transform the equation, and then work with the transformed equation. In effect, one has to write a new code to analyze the old problem. The original reason we were attracted to the template approach to factoring out translational invariance was that we were interested in using it in what in general we term ”time-stepper” based methods for numerical bifurcation theory. It is possible to implement the template procedure not only in continuous time, but also in discrete time: one runs the “original” dynamics for some time and then “pulls back” the result (as opposed to constantly pulling it back as in ). We have implemented this discrete time, time-stepper based bifurcation approach in  for the traveling case; in the same spirit, the discrete-time version of the approach can be used on the original equation in the self-similar case (as opposed to writing down the transformed MN-equation).

This opens the way for another interesting possibility. Over the last few years, we have developed methods for what we call “coarse bifurcation analysis” of microscopic time-steppers . In these, we create a map (using microscopic dynamics, such as Monte Carlo, or Lattice Boltzmann) from an initial condition in the space of a few moments of a distribution to a final condition (after some time) in the same truncated moments space of the distribution. The procedure involves “lifting” the macroscopic initial condition to one or more microscopic ones conditioned on the few governing moments, running the microscopic code, and then averaging (restricting) back to governing moments space. There is a number of interesting issues about how long to run the microscopic simulator, and how to do good variance reduction, but we will not discuss this here.

Based now on the coarse-time stepper implementation, suppose that we have microscopic simulations (simulations at one level of description), and we suspect that some coarser description (moments) of the microscopic dynamics will have self-similar behavior. Then, using the coarse time-stepper and the discrete-time version of the above procedure, it is possible, given a set of reasonable assumptions discussed in some detail in  to find certain elements of the self-similar dynamics (e.g. steady states, invariant sets, bifurcations) of what the “coarse-MN-equation” would have been, without having to explicitly construct (approximate) this coarse-MN equation. Indeed, this approach gives rise to alternative, mathematics-motivated ensembles for performing the microscopic simulations (with an eye towards the analysis of coarse dynamics) . Recent advances in multiscale computations (projective and telescopic projective integrators , as well as what we provisionally call “patch dynamics”) may assist in doing these calculations efficiently, in smaller space and time scales than the full computational domain. These developments in multiscale computations are rather general, and not designed in particular with the coarse MN-dynamics in mind.

We believe that the general approach outlined in this paper opens new computational possibilities in the study of self-similar problems and their dependence on parameters. Problems ranging from interfacial fluid dynamics to shock waves, and from cell aggregation to cosmology and materials science may become more accessible to computer-assisted analysis, whether through closed macroscopic equations, or through the “coarse” analysis of alternative, more microscopic descriptions.

## 6 Acknowledgements

I.G.K. would like to acknowledge partial support by the AFOSR (Dynamics and Control, Dr. Jacobs and Dr. King); also the visit of Clancy Rowley (then a Caltech graduate student, now a colleague) to MAE/PACM in Princeton in 1999, and the collaboration with Drs. Runborg and Theodoropoulos in implementing the Rowley-Marsden approach for discrete time “effective equation” problems .


text:
# 1 Introduction

## 1 Introduction

Violation of CP, T and CPT symmetries in the time evolution of $`K^0`$-$`\overline{K}^0`$ was studied by the CP-LEAR experiment  from the preparation of definite flavour states. The study of this *flavour-to-flavour* evolution allows the construction of observables which violate CP and T, or CP and CPT. In order to be non-vanishing, nevertheless, these observables need the presence of an absorptive part in the effective Hamiltonian that governs neutral meson system. The different lifetimes of physical states $`K_L`$ and $`K_S`$ provides this ingredient. In the case of $`B_d`$ mesons, on the contrary, the width difference $`\mathrm{\Delta }\mathrm{\Gamma }`$ between the physical states is expected to be negligible, so that the T- and CPT-odd observables proposed for kaons, and based on flavour tag, will practically vanish for a $`B_d`$ system.

Here alternative observables are discussed, which allow the study of CP, T and CPT indirect violation in the $`B_d`$ system. Based on CP-tag, these observables do not need the presence of $`\mathrm{\Delta }\mathrm{\Gamma }0`$, and can be constructed from the entangled states of $`B_d`$ mesons.

In the following section, the invariant parameters $`\epsilon `$ and $`\delta `$ are introduced to describe indirect violation of symmetries in the neutral meson system. In section 3 we describe the CP tag of $`B_d`$ from the entangled states in a $`B`$-factory. Next, section 4 reviews three different kinds of asymmetries that can be constructed from these states, namely, *flavour-to-flavour* and both genuine and non-genuine *CP-to-flavour* asymmetries. Finally, in section 5 the first estimates on the reach and sensitivity of the experimental analysis are given.

## 2 Invariant description of CP, T and CPT violation in the $`B`$ system

The physical states in the neutral $`B`$-meson system are a linear combination of the definite flavour $`B^0`$ and $`\overline{B}^0`$. Physical states can also be written in terms of CP eigenstates, $`|B_\pm \frac{1}{\sqrt{2}}(I\pm CP)|B^0`$, which are physical iff the CP operator is well defined. To do so, one has to introduce two complex parameters, $`\epsilon _{1,2}`$, to describe the CP mixing, so that $`|B_{1(2)}=\frac{1}{\sqrt{1+|\epsilon _{1(2)}|^2}}\left[|B_{+()}+\epsilon _{1(2)}|B_{(+)}\right]`$. The complex parameters $`\epsilon _{1,2}`$, invariant under rephasing of the meson states, are better interpreted in terms of $`\epsilon (\epsilon _1+\epsilon _2)/2`$ and $`\delta \epsilon _1\epsilon _2`$, whose observable character is explicit when they are written in terms of the effective hamiltonian matrix elements .

Discrete symmetries impose different restrictions on the effective mass matrix, $`H=M\frac{i}{2}\mathrm{\Gamma }`$, and thus on the invariant parameters $`\epsilon `$ and $`\delta `$:

* CPT invariance requires222Here $`H_{ij}`$, $`M_{ij}`$, and so on, represent the matrix elements in the flavour basis $`B^0\overline{B}^0`$. $`H_{11}=H_{22}`$, so that $`\delta =0`$, with no restriction on $`\epsilon `$;
* T invariance imposes $`\mathrm{Im}(M_{12}\mathrm{CP}_{12}^{})=\mathrm{Im}(\mathrm{\Gamma }_{12}\mathrm{CP}_{12}^{})=0`$, and so $`\epsilon =0`$;
* and CP conservation requires both $`\epsilon =\delta =0`$.

In the exact limit $`\mathrm{\Delta }\mathrm{\Gamma }=0`$, an approximation that is expected to be excellent for the $`B_d`$ system, both $`\mathrm{Re}(\epsilon )`$ and $`\mathrm{Im}(\delta )`$ vanish. Then $`\mathrm{Im}(\epsilon )0`$ is a proof of both CP and T violation, and $`\mathrm{Re}(\delta )0`$ is a proof of CP and CPT violation, but neither $`\mathrm{Re}(\epsilon )=0`$ nor $`\mathrm{Im}(\delta )=0`$ are proof of a fundamental invariance. Information on the symmetry parameters can be extracted from the study of time evolution of $`B`$ meson entangled states.

## 3 CP-Tag from entangled states

In a $`B`$ factory operating at the $`\mathrm{{\rm Y}}(4S)`$ peak, correlated pairs of neutral $`B`$-mesons are produced through $`e^+e^{}\mathrm{{\rm Y}}(4S)B\overline{B}`$. Charge conjugation together with Bose statistics require the initial state to be

$$|i>=\frac{1}{\sqrt{2}}\left(|B^0(\stackrel{}{k}),\overline{B}^0(\stackrel{}{k})>|\overline{B}^0(\stackrel{}{k}),B^0(\stackrel{}{k})>\right).$$
(1)
This permits the performance of a flavour tag: if at $`t_0`$ one of the mesons decays through a channel $`X`$, which is only allowed for one flavour, the other meson in the pair must have the opposite flavour at $`t_0`$, and will later evolve during $`\mathrm{\Delta }t=tt_0`$ until its final decay to some state $`Y`$.

The entangled $`B\overline{B}`$ state can also be expressed in terms of the CP eigenstates as $`|i>=\frac{1}{\sqrt{2}}\left(|B_{},B_+>|B_+,B_{}>\right)`$. Thus it is also possible to carry out a CP tag, if we have a CP-conserving decay into a definite CP final state $`X`$, so that its detection allows us to identify the decaying meson as a $`B_+`$ or a $`B_{}`$, which decays into $`Y`$ after a time $`\mathrm{\Delta }t`$. In Ref.  we described how this determination is possible and unambiguous to $`𝒪(\lambda ^3)`$, the flavour-mixing parameter of the CKM matrix.

If we consider only decay channels $`X`$, $`Y`$ which are either flavour or CP conserving, then the final configuration $`(X,Y)`$ corresponds to a single particle mesonic transition. The intensity for the final configuration, $`I(X,Y,\mathrm{\Delta }t)\frac{1}{2}_{\mathrm{\Delta }t}^{\mathrm{}}𝑑t^{}|(X,Y)|^2`$ is proportional to the time dependent probability for the meson transition.

## 4 Asymmetries

By comparing the probabilities corresponding to different processes we build time-dependent asymmetries that can be classified into three types.

### 4.1 Flavour-to-flavour genuine asymmetries

The final configuration denoted by $`(\mathrm{},\mathrm{})`$, with flavour definite (for example, semileptonic) decays detected on both sides of the detector, corresponds to *flavour-to-flavour* transition at the meson level. The equivalence is shown in Table 1.

The first two processes in the Table are conjugated under CP and also under T. The corresponding Kabir asymmetry is, to linear order in the CPT violating $`\delta `$,

$$A(\mathrm{}^+,\mathrm{}^+){\scriptscriptstyle \frac{4\frac{\mathrm{Re}(\epsilon )}{1+|\epsilon |^2}}{1+4\frac{\mathrm{Re}(\epsilon )}{1+|\epsilon |^2}}},$$
(2)
which does not depend on time. However, in the exact limit $`\mathrm{\Delta }\mathrm{\Gamma }=0`$, $`\mathrm{Re}(\epsilon )`$ vanishes, and this quantity will be zero. For the $`B_d`$ system, experimental limits on $`\mathrm{Re}(\epsilon )`$ are of few parts in a thousand .

A second asymmetry arises from the last two processes in Table 1, related by a CP or a CPT transformation,

$$A(\mathrm{}^+,\mathrm{}^{})2\frac{\mathrm{Re}({\scriptscriptstyle \frac{\delta }{1\epsilon ^2}})\mathrm{sinh}{\scriptscriptstyle \frac{\mathrm{\Delta }\mathrm{\Gamma }\mathrm{\Delta }t}{2}}\mathrm{Im}({\scriptscriptstyle \frac{\delta }{1\epsilon ^2}})\mathrm{sin}(\mathrm{\Delta }m\mathrm{\Delta }t)}{\mathrm{cosh}{\scriptscriptstyle \frac{\mathrm{\Delta }\mathrm{\Gamma }\mathrm{\Delta }t}{2}}+\mathrm{cos}(\mathrm{\Delta }m\mathrm{\Delta }t)},$$
(3)
which is an odd function of time. This asymmetry also vanishes unless $`\mathrm{\Delta }\mathrm{\Gamma }0`$. Present limits  on $`\mathrm{Im}(\delta )`$ are at the level of few percent.

### 4.2 CP-to-flavour genuine asymmetries

Alternative asymmetries can be constructed making use of the CP eigenstates, which can be identified in this system by means of a CP tag. If the first decay product, $`X`$, is a CP eigenstate produced along the CP-conserving direction, i.e. the decay is free of CP violation, and $`Y`$ is a flavour definite channel, then the mesonic transition corresponding to the configuration $`(X,Y)`$ is of the type *CP-to-flavour*.

In Table 2 we show the mesonic transitions, with their related final configurations, connected by genuine symmetry transformations to $`B_+B^0`$.

Comparing the intensities of the four processes, we may construct three genuine asymmetries, namely $`A(\mathrm{CP})`$, $`A(\mathrm{T})`$ and $`A(\mathrm{CPT})`$ .

$$A(\mathrm{CP})=2{\scriptscriptstyle \frac{\mathrm{Im}(\epsilon )}{1+|\epsilon |^2}}\mathrm{sin}(\mathrm{\Delta }m\mathrm{\Delta }t)+{\scriptscriptstyle \frac{1|\epsilon |^2}{1+|\epsilon |^2}}{\scriptscriptstyle \frac{2\mathrm{R}\mathrm{e}(\delta )}{1+|\epsilon |^2}}\mathrm{sin}^2\left({\scriptscriptstyle \frac{\mathrm{\Delta }m\mathrm{\Delta }t}{2}}\right),$$
(4)
the CP odd asymmetry, contains both T-violating and CPT-violating contributions, which are, respectively, odd and even functions of $`\mathrm{\Delta }t`$. This asymmetry corresponds to the ”gold plate” decay  and has been measured recently . T and CPT violating terms can be separated by constructing other asymmetries.

$$A(\mathrm{T})=2{\scriptscriptstyle \frac{\mathrm{Im}(\epsilon )}{1+|\epsilon |^2}}\mathrm{sin}(\mathrm{\Delta }m\mathrm{\Delta }t)\left[1{\scriptscriptstyle \frac{1|\epsilon |^2}{1+|\epsilon |^2}}{\scriptscriptstyle \frac{2\mathrm{R}\mathrm{e}(\delta )}{1+|\epsilon |^2}}\mathrm{sin}^2\left({\scriptscriptstyle \frac{\mathrm{\Delta }m\mathrm{\Delta }t}{2}}\right)\right],$$
(5)
the T asymmetry, needs $`\epsilon 0`$, and turns out to be purely odd in $`\mathrm{\Delta }t`$ in the limit we are considering.

$$A(\mathrm{CPT})={\scriptscriptstyle \frac{1|\epsilon |^2}{1+|\epsilon |^2}}{\scriptscriptstyle \frac{2\mathrm{R}\mathrm{e}(\delta )}{1+|\epsilon |^2}}\frac{\mathrm{sin}^2\left({\scriptscriptstyle \frac{\mathrm{\Delta }m\mathrm{\Delta }t}{2}}\right)}{12\frac{\mathrm{Im}(\epsilon )}{1+|\epsilon |^2}\mathrm{sin}(\mathrm{\Delta }m\mathrm{\Delta }t)},$$
(6)
is the CPT asymmetry. It needs $`\delta 0`$, and includes both even and odd time dependences.

The above expressions correspond to the limit $`\mathrm{\Delta }\mathrm{\Gamma }=0`$, but, being genuine observables, a possible absorptive part could not induce by itself a non-vanishing asymmetry.

### 4.3 CP-to-flavour non-genuine asymmetries

The construction of the quantities described in the previous paragraphs requires to tag both $`B_+`$ and $`B_{}`$ states, and thus the reconstruction of the experimentally challenging decay $`BJ/\mathrm{\Psi }K_L`$. Conversely, non-genuine asymmetries offer a possibility to measure the symmetry parameters from the reconstruction of $`J/\mathrm{\Psi }K_S`$ only. But they involve the discrete transformation that we denote $`\mathrm{\Delta }t`$, consisting of the exchange in the order of appearance of decay products $`X`$ and $`Y`$, which cannot be associated with any fundamental symmetry.

Table 3 shows the different transitions we may study from such final states.

Besides the genuine CP asymmetry, there are two new quantities that can be constructed from the comparison between $`(J/\mathrm{\Psi }K_S,\mathrm{}^+)`$ and the processes in the table. In the exact limit $`\mathrm{\Delta }\mathrm{\Gamma }=0`$, $`\mathrm{\Delta }t`$ and T operations are found to become equivalent, so that the temporal asymmetry satisfies $`A(\mathrm{\Delta }t)A(\mathrm{}^+,J/\mathrm{\Psi }K_S)=A(\mathrm{T})`$ and moreover $`A(\mathrm{CP}\mathrm{\Delta }t)A(\mathrm{}^{},J/\mathrm{\Psi }K_S)=A(\mathrm{CPT})`$. Since this result holds for $`\mathrm{\Delta }\mathrm{\Gamma }0`$, it is expected to be valid for the $`B_d`$ system, but not for $`B_s`$ and even less for $`K`$. The asymmetries $`A(\mathrm{\Delta }t)`$ and $`A(\mathrm{CP}\mathrm{\Delta }t)`$ are non-genuine, and the presence of $`\mathrm{\Delta }\mathrm{\Gamma }0`$ may induce non-vanishing values for them, even if there is no true T or CPT violation. These fake effects, nevertheless, can be calculated and are thus controllable.

## 5 CP, T, CPT indirect violation reach at asymmetric B-Factories

The asymmetries described in the previous section can be already constructed from the current data taken at Asymmetric B-Factories . The experimental analysis is based on a simultaneous unbinned likelihood fit of the flavour and CP intensities $`I(X,Y;\mathrm{\Delta }t)`$, together with the $`B^0`$/$`\overline{B}^0`$ mistag rates and the $`\mathrm{\Delta }t`$ resolution function. The coefficients of terms with different temporal dependencies contain the information on the symmetry parameters.

From a simulation study, estimations on the reachable statistical precision for the relevant parameters have been calculated for $``$60 fb-1 (assuming yields from Ref. ) and are shown in Table 4.

## 6 Conclusions

We have shown how the two complex rephasing invariant parameters $`\epsilon `$ and $`\delta `$ describe CP, T and CPT indirect violation in $`B^0\overline{B}^0`$. In the exact limit $`\mathrm{\Delta }\mathrm{\Gamma }=0`$ the number of parameters is reduced to $`\mathrm{Im}(\epsilon )`$ and $`\mathrm{Re}(\delta )`$. Observables based on flavour-to-flavour transitions are sensitive to $`\mathrm{Re}(\epsilon )`$, but need $`\mathrm{\Delta }\mathrm{\Gamma }0`$, and thus are not promising in B-factories. Conversely, these experimental facilities allow the construction of new asymmetries based on combination of flavour and CP tags.

First estimations on the sensitivity reachable on B-factories have been presented. This data will be crucial to achieve the separation of the two ingredients: on one hand CP and T violation, described by $`\epsilon `$, and on the other CP and CPT violation, given by $`\delta `$.

This work has been supported by CICYT, Spain, under Grant AEN99-0692.


text:
# The role of chaotic resonances in the solar system

Our understanding of the solar system has been revolutionized by the finding that multiple planet systems are subject to chaotic dynamical processes. The giant planets in our solar systems are chaotic, as are the inner planets (independently). Recently a number of extrasolar planetary systems have been discovered; some of these systems contain multiple planets. It is likely the some of the latter are chaotic. In extreme cases chaos can disrupt some orbital configurations, resulting in the loss of a planet. The spin axes of planets may also evolve chaotically. This may adversely affect the climate of otherwise biologically interesting planets.

Observers have targeted solar-type stars in their searches for planets, partly motivated by the belief that these stars might harbor Earth-like planets. Although present detection techniques are unable to detect terrestrial-mass planets, they are sensitive to Jupiter-mass bodies. The original notion was that because all of these solar systems should share a common planet formation process they would result in similar collections of planets; large planets might be accompanied by smaller counterparts. However, the wide variety of orbits of the planets detected to date shows that our knowledge of planet formation is incomplete. Nevertheless, while the formation processes in other solar systems might be quite different from those that shaped our own, the long-term dynamical evolution is governed by the same principles.

These principles have been known since Newton, but their consequences for our solar system still surprise researchers. For nearly 300 years astronomers believed that the orbits of the planets were regular and predictable. After an initial phase of formation, the number of planets, asteroids, and comets was fixed. Planets neither escaped the solar system nor collided with each other. The discovery of chaos destroyed this traditional picture.

Our solar system provides a plethora of examples of chaotic motion. The theory of chaos has been used to explain, and in some cases predict, the location and extent of gaps in the asteroid and Kuiper belts1 -7 . The theory predicts that irregularly-shaped satellites such as Hyperion tumble chaotically8 . The obliquity of Mars undergoes large excursions9 ; these are in fact chaotic10 ,11 . Even the obliquity of the distant future Earth may undergo chaotic evolution in 1.5-4.5 Gyr12 -14 . Furthermore, we have recently shown that the chaos in the orbits of the giant planets, first observed by Sussman and Wisdom15  is due to a three-body resonance among Jupiter, Saturn, and Uranus16 ; the locations of these planets cannot be predicted on a time scale longer than a few tens of millions of years. This provides analytic support to earlier numerical experiments that demonstrated that our planetary system evolves chaotically15 ,17 . The worried reader may find some comfort; the accompanying analytic theory predicts that no planet will be ejected before the sun dies.

The rigid pendulum

Despite the variety and complexity of the applications mentioned above, we can introduce many of the concepts in solar system dynamics using the pendulum: phase space structure, periodic motion, and stability. To describe the state of the pendulum we must specify both its position (such as the angle from the downward vertical, $`\theta `$, see figure 1) and velocity (such as the time rate of change of the angle, $`\dot{\theta }`$, or the momentum $`p=ml\dot{\theta }`$, where $`m`$ is the effective mass of the pendulum and $`l`$ is the length). The structure of the phase space is built around fixed points (see figure 2). When the pendulum is hanging straight down and motionless we have a fixed point. There is another fixed point when the pendulum is straight up and motionless. Readers equipped with a pencil can immediately see that these two points are different: dangling a pencil loosely between two fingers is easy, but balancing a pencil upright challenges the most dextrous.

If a pendulum at the downward fixed point is nudged it begins to oscillate or librate back and forth; this fixed point is said to be stable. If a pendulum momentarily balanced on its head is gently nudged it swings down and back up on the other side; this fixed point is unstable. If the pendulum is pushed forcefully enough it will rotate periodically. Thus, we see three main regions of the phase space: circulation in one direction, libration, and circulation in the other direction. These three regions are separated by the curve that connects the unstable fixed point to itself. This curve forms the “cat’s eye” in figure 2, and is called the separatrix. All trajectories inside the separatrix are said to be in resonance. Trajectories well outside the “cat’s eye” have enough kinetic energy that they are unaffected by the gravitational potential. They rotate with an almost fixed angular frequency. In contrast trajectories inside the resonance are dominated by the potential; they librate. The average value of $`\dot{\theta }`$ in the resonance is zero, which is a generic signature of resonant phenomena.

The separatrix is the cradle of chaos. Suppose we reach in a hand and just nudge a moving pendulum in an attempt to change the character of its motion. If the pendulum is far from the separatrix (either librating near the stable fixed point or rotating) the nudge will not change the character of the motion. However, close to the separatrix the pendulum nearly stands motionless on its head once every libration or rotation period. A nudge can then change rotation into libration, libration into rotation, or rotation in one direction into rotation in the other direction. Aside from the separatrix, all the orbits in figure 2 lie on closed curves. These curves are called invariant curves; similar curves appear in most dynamical systems.

In order for chaos to occur there must be at least two interacting oscillators. In the solar system that interference is supplied by a third body, the invisible hand made immanent. This interference manifests itself as zones of chaotic motion, which are the ghosts of the departed separatrices. That phase space would harbor regions where trajectories are largely quasiperiodic and regions where trajectories are largely chaotic is a feature of generic dynamical systems. This characteristic is known as a divided phase space18 ,19 .

Kirkwood gaps

A series of remarkable features in the asteroid belt vividly illustrates the importance of dynamical chaos in the solar system. The distribution of semimajor axes of the orbits of the asteroids contains a number of distinct gaps. These are call Kirkwood gaps, in honor of Daniel Kirkwood20 . It was he who first identified them and noted that they occur at locations where the orbital period, $`T`$, which depends on semimajor axis, would be of the form $`(p/q)T_J`$, where $`T_J`$ is the orbital period of Jupiter and $`p`$ and $`q`$ are integers. In terms of the orbital frequencies, $`n=2\pi /T`$, $`qn_Jpn0`$; in other words this is a resonance. Nearly one hundred years passed it was explained how a resonance produced a shortage of asteroids, or a gap. The reason for this delay can be appreciated by making a simple estimate of the strength of a resonant perturbation. The force is proportional to the mass of Jupiter, which is $`\mu _J0.001`$ in units of the solar mass. Worse still, the resonant force is zero if the orbits are circular (for a planar model). If the orbits are non-circular, with an eccentricity $`e`$ for the asteroid and $`e_J`$ for Jupiter, the force is proportional to $`e^{|pq|s}e_J^s`$, where $`s`$ is an integer between zero and $`|pq|`$. For example, for a 3:1 resonance $`p=3`$ and $`q=1`$, three force terms are possible: $`e^2`$, $`ee_J`$, and $`e_J^2`$. Since $`ee_J0.05`$, the force exerted on the asteroid by the resonance is smaller by a factor of $`5\times 10^5`$ than the force exerted by the sun, in the best possible case of a first order resonance ($`|pq|=1`$).

Laplace noted even before the discovery of the Kirkwood gaps that a resonant force adds up over many orbits, but only up to half the libration period (the time to go half way around the cat’s eye in figure 2) which is proportional to the inverse of the square root of the resonant force. The integrated force is then larger, roughly like the square root of $`\mu _Je^{|pq|}`$. This is still a small effect, so it is difficult to see how the asteroid can be ejected from the resonance.

Even if the force added up over longer times, theorems on the behavior of dynamical systems from the 1950’s and ’60’s suggested that gap formation might be difficult. For example, the motion of a typical asteroid in an imaginary solar system with very small planetary masses lies on invariant curves similar to those in figure 2 (see Box 1). This type of result is known as a KAM theorem21 -23 . In exceptional cases, near the separatrix of a resonance, the theorem does not apply.

The fact that the theorems are restricted to very small masses and the demonstration of a divided phase space suggest a way to evade the dilemma; perhaps for realistic values, invariant curves at substantial distances from the separatrix are destroyed. Numerical integrations by Giffen of a simplified model, incorporating only a single planet (Jupiter) on a fixed ellipse, showed that the invariant curves were destroyed in the vicinity of the 2:1 resonance1 . This removed the barrier of the KAM theorem. Unfortunately, the $``$100,000 yr integrations by Giffen and others, while establishing the presence of chaotic orbits, did not give any hint as to how the irregular motion led to the removal of an asteroid from the resonance1 ,24 .

The critical advance came when Wisdom showed by direct numerical integration that the eccentricity of small bodies placed in the 3:1 gap alternates chaotically between periods of low and moderate eccentricity as the result of perturbations from Jupiter2  (see figure 4). This result relied upon upon the development of an algebraic mapping that could efficiently follow the long-term motion of asteroids in the 3:1 mean motion resonance with Jupiter for a million years or longer. This is a recurring theme in solar system dynamics: more efficient algorithms and faster computers permit longer, more accurate integrations often reveal unexpected dynamical results.

What is the source of the chaos in the 2:1 and 3:1 mean motion resonances? It has long been known that each mean motion resonance is composed of several individual resonant terms25 . For example in the 3:1 there are three such terms in the planar case, as mentioned above. While each of these resonances are weak, they are also close together2 -4 ,26 . It is the close proximity of multiple, albeit weak, resonances which invalidates the simple estimate for the integrated force given above. A single resonance produces a small force integrated over a libration period, and a simple “cats eye” in phase space with a smooth separatrix. Adding a second weak resonance breaks up this smooth separatrix in a chaotic region (see Box 1). The separatrix is broken because the second resonance nudges the pendulum represented by the first resonance, altering the motion, particularly of orbits that pass near the unstable fixed point50 . As a result the motion from one passage near the fixed point to the next is uncorrelated. The resonant force is still small compared to the force exerted by the sun, and produces only small changes over a libration period. But over times longer than a libration period the asteroid experiences uncorrelated forces.

As a result of these uncorrelated forces the motion in the chaotic zone differs from that in a regular region; the chaotic motion is very sensitive to initial conditions. Suppose we take two copies of our interacting oscillators, released from slightly different positions. In one case we start near a stable fixed point, well away from the separatrix. Here the motion is regular and the two copies of the apparatus separate from each other linearly (or polynomially) in time (see figure 5). In the other case we start the two copies near an unstable fixed point. Now the separation increases exponentially in time. The time scale for this separation is called the Lyapunov time, $`T_L`$.

We can estimate the Lyapunov time of an asteroid placed artificially in the 3:1 mean motion resonance. There are only two relevant time scales in the problem: the libration period and the precession period. From our argument above, large-scale chaos occurs when the separation of the islands is comparable to their widths27 , that is, when the resonances overlap. Another way to say this is that large scale chaos occurs when the libration and precession periods of a resonant asteroid are similar. In that situation the Lyapunov time must be equal to the libration period—the only time scale in the problem28 . In the case of the 3:1 resonance the precession period is $`3\times 10^4`$ years while the numerically-determined Lyapunov time is $`1.4\times 10^4`$ years.

The two types of trajectories differ in more than their rates of separation. The chaotic trajectory can explore a larger volume of phase space than is accessible by a regular trajectory. As noted by Wisdom, for asteroids this extra measure of freedom can be dangerous, opening routes by which they can be ejected or cannibalized by planets2 -5  (it’s world-eat-world out there).

So far we have discussed regular and chaotic motion as if the two occurred in different systems, but in fact most dynamical systems exhibit a rich mixture of regular and chaotic motion. Box 1 describes this in more detail.

In higher order resonances, where the phase space structure is relatively simple, the uncorrelated or chaotic forces lead to simple random walks involving only small exchanges of energy between the planet and the asteroid, since the semimajor axes of both bodies are nearly fixed by the resonance condition. Not so for the angular momentum; the torques associated with the uncorrelated forces produce a random walk in the angular momentum of both bodies. This leads to a classic gambler’s ruin problem with the asteroid playing the part of the gambler and the planet taking the role of the house. If either body loses most of its angular momentum, that body’s orbit becomes highly eccentric and subject to collisions with other objects in the solar system or ejection—the analog of bankruptcy. The planet, with its much larger mass, has more capital (angular momentum) than the asteroid. Even worse, the asteroid is like a gambler forced to limit his winnings to what he can carry in his pockets; because of its small mass the asteroid cannot absorb enough angular momentum to produce a substantial change in the orbit of the planet. The two trade angular momentum back and forth, but when the inevitable happens and the asteroid loses its small stock of angular momentum, it is abruptly escorted from the cosmic casino.

We can estimate the time to remove the asteroid. Because the motion is chaotic we can treat the resonance angle as a random variable, which drives random changes in the eccentricity. In the case of the 3:1 resonance the changes in $`e`$ are of order $`\delta e^2(T_L/T_J)(M_J/M_{})e^2`$, using the Lyapunov time, $`T_L`$, as the interval associated with the random changes ($`M_J`$, $`T_J`$, and $`M_{}`$ are Jupiter’s mass, Jupiter’s orbital period, and the solar mass, respectively). The short-term average eccentricity diffuses to large values in a time of about one million years.

In some cases, such as the asteroid Helga in the 12:7 resonance, the motion is chaotic with a very short Lyapunov time, but the diffusion time is comparable to or larger than the age of the solar system. When applied to Helga this theory predicts that the asteroid should survive for something like 8 billion years. Direct numerical integrations agree with this prediction6 . We look forward to seeing this prediction verified observationally.

It is worth noting that the largest gap, that between the 1:1 resonance harboring the Trojan asteroids, and the Hildas in the 3:2 resonance, is the result of the overlap of distinct mean motion resonances29 ; for example the 4:3 resonance overlaps with the 5:4, the 5:4 with the 6:5, and so forth.

It is believed that most meteorites come from one of two sources, the 3:1 mean motion resonance and the $`\nu _6`$ secular resonance4 ,30 ,31 . The latter involves a resonance between the precession frequency of the apsidal line of the asteroid and the sixth fundamental secular frequency of the solar system (which is very roughly Saturn’s precession frequency). Roughly equal numbers of meteorites come from each type of resonance. Until recently there was one difficulty with this story. The cosmic ray exposure ages (a measure of the delivery time) of the stony meteorites are typically 20 Myr32 , where as the delivery time from the 3:1 resonance is much shorter, 1 Myr, as noted above. However, recent simulations by Vokruhlický and Farinella33  suggests a way out of this dilemma. They suggest that most meteorites, which are produced by collision between larger bodies in the asteroid belt, are not injected directly into either the 3:1 or the $`\nu _6`$ resonance. Rather they are placed in the vicinity of the resonance, and then dragged into the unstable region by the Yarkovsky effect, which arises due to the anisotropic thermal radiation from those regions of the fragment which are exposed to the sun. With the addition of this slow dynamical precursor, the story of meteorite delivery appears to be complete.

Chaotic Spin-Orbit Resonances

One of the most dramatic and important examples of chaos is afforded by the evolution of planetary spins. This chaos is produced by resonances between spin and orbit precession periods. In such resonances, the asphericity of the planet couples to the non-axisymmetric perturbation produced by orbital eccentricity or inclination. The rather small deviations from perfect sphericity, typically a few parts in a thousand for larger bodies like planets, leads to significant exchanges of energy and angular momentum between orbital and rotational (or spin) motions of satellites and even planets. Resonances between spin period and orbital period are common in the solar system, the Moon being a prominent example. The phase space occupied by such resonances is small, making it unlikely that so many bodies formed in resonance. This paradox is explained by the fact that dissipative effects tend to drive bodies into these resonances, where the motion is stablized34 ,35 . As a result, most resonant satellites are currently deep in their respective 1/1 spin-orbit resonances, where the motion is regular; however, in one case, Saturn’s satellite Hyperion, the chaotic motion may have persisted until the present8 .

Because solar tides are so weak, dissipative effects tend to be less important for planetary spins, allowing for richer dynamics. Seminal work by Ward36  showed that the angle between the spin and orbital axes of Mars (the “obliquity”) varies by a $`\pm 13.6^{}`$ around its average of $`24^{}`$ over millions of years. These variations are a result of a resonance between the precession of the spin axes and a combination of orbital precession frequencies; improvements in orbital models resulted in an increased variation of $`\pm 20^{}`$9 . A dramatic twist was added by numerical integrations, which showed that the obliquity of Mars is evolving chaotically, and varies over an even larger range10 ,11 . Such a variation has profound, but as yet poorly understood implications for climate variation. It is likely that the spin axes of Mercury and Venus underwent chaotic variations in the past11 .

The tilt of the Earth, currently 23 degrees, will also increase in the future, as first pointed out by Ward12 ; the Moon evolves outward under the influence of the tides resulting in a decrease in the precession rate of the Earth. Eventually the precession rate becomes resonant with yet another combination of orbital precession rates. Once again, numerical integrations show that the Earth’s obliquity will vary chaotically13 ,14 . The tilt of the Earth’s axis may increase to 90 degrees. The effect on our climate is hard to estimate, but the result is unlikely to be pleasant.

Three-Body Resonances

Numerical integrations of main belt asteroids show that a substantial fraction of these bodies have chaotic orbits with rather short Lyapunov times ($`10^5`$ years)6 ,37 ,38 . This chaotic motion is not associated with any two-body resonance. Instead it is the result of the interaction among three bodies: the asteroid, Jupiter, and either Mars in the inner asteroid belt or Saturn in the main or outer asteroid belt37 ,38 . As with Helga, the Lyapunov time can be quite short, but the diffusion times are comparable to or longer than the age of the solar system.

Three-body resonances arise when one planet, Jupiter say, is perturbed by a second, Saturn for example. The orbit of Jupiter is then no longer a simple Keplerian ellipse. The potential experienced by the asteroid in Jupiter’s gravitational field is given by an expression formally equivalent to the two-body case, but Jupiter’s orbital elements now vary with time. This variation introduces a whole new suite of frequencies into the potential experienced by the asteroid; in addition to all the harmonics of Jupiter’s period, all the harmonics of Saturn’s orbital period appear as well. Three-body resonances have also been considered for the Galilean and Uranian satellites, as well as in ring systems39 -41 .

These new terms in the potential have much smaller amplitudes than two body resonant terms having the same number of powers of eccentricity. Three-body resonances are proportional to the product of the masses of the two perturbing bodies, $`M_J`$ and $`M_S`$ in our example, rather than just $`M_J`$ in a two body resonance.

The structure of a three-body resonance is similar to that of a two body resonance—multiple, very narrow components separated in semimajor axis by an amount that depends on the precession frequencies of the bodies involved. When the separation is comparable to or smaller than the width of the individual resonances, the motion in the immediate vicinity of the resonance will be chaotic. One can then estimate the Lyapunov time and the diffusion time in the same manner as for two body resonances.

We have searched the catalogue of asteroids for evidence of gaps at the location of a number of the stronger three-body resonances. We have not been able to identify a candidate gap.

We can check the theory in other ways. For example, some asteroids on the inner edge of the belt (near Mars) could have formed in three-body resonances involving the asteroid, Mars, and Jupiter. If the diffusion time of some of these asteroids is comparable to the age of the solar system, we should be able to find objects that are about to be removed from the belt due to close encounters, or possibly collisions, with either Mars or Earth. Morbidelli and co-workers claim that such bodies have already been discovered, namely the so-called near-Earth asteroids42 . These bodies have been a puzzle for astronomers, since in their present orbits they have very short lifetimes, typically millions of years. Since the solar system is billions of years old, objects with lifetimes of millions of years should have vanished long ago, unless there is some way to replenish the supply. The theory of three-body resonances offers a possible mechanism. In contrast to the situation with Helga, we would prefer not to test the predictions of the theory, at least in its collisional aspects.

Chaos among the giant planets

Even trajectories outside but near a resonance can be affected by its presence. In the century following Newton’s publication of his law of universal gravitation, astronomers noted that the positions of Jupiter and Saturn deviated from their predicted positions by some 30 minutes of arc. The difference became known as the great inequality. Laplace noted that the predictions did not take into account the influence of Saturn on Jupiter’s orbit, and vice versa. On the face of it this neglect seemed reasonable, because the mass of either planet was less than one one-thousandth than that of the sun. Early astronomers had employed the naive estimate made above, which integrates this tiny force over the interval between successive conjunctions of the planets, a time somewhat longer than $`T_J`$. But Laplace realized that the orbital period of Saturn was almost exactly $`5/2`$ times that of Jupiter ($`\frac{2}{T_J}\frac{5}{T_S}\frac{1}{85T_J}`$). The perturbations accumulate over a much longer interval ($`85T_J`$), permitting a larger exchange of energy and angular momentum between the two planets. The change in the predicted position of the planet on the sky was roughly $`85^2`$ times larger than the simple estimate would indicate. With this result in hand, Laplace was able to reconcile the observations with the prediction of the law of gravity.

This discovery strongly affected Laplace’s views regarding determinism, reflected in his well known statement43 :

> The present state of the system of nature is evidently a consequence of what it was in the preceding moment, and if we conceive of an intelligence that at a given instant comprehends all the relations of the entities of this universe, it could state the respective position, motions, and general affects of all these entities at any time in the past or future.
>
> Physical astronomy, the branch of knowledge that does the greatest honor to the human mind, gives us an idea, albeit imperfect, of what such an intelligence would be.

This view of the world passed into common currency: the clockwork motion of the planets became the epitome of regularity.

This view is wrong. It is wrong because it ignores the delicate nature of the separatrix, the cradle of chaos. In fact the giant planets provide further evidence of the power of three-body resonances, the subtlety of natural phenomena, and the difficulty of interpreting the fruits of scientific enquiry.

The tour-de-force numerical integrations of the outer planets by Sussman and Wisdom in 198844  shattered the clockwork. Using the “Digital Orrery”, a parallel computer built specifically for the the task and now part of the Smithsonian collection, they followed the motion of the four giant planets and Pluto (as a test particle) over 845 Myr. To the surprise of all at the time, Pluto’s orbit was chaotic with $`T_L10`$ million years.

This revolutionary breakthrough was followed by a re-examination of the orbital evolution of the full solar system. Longer integrations 1989 by Laskar of an approximate model of the solar system (the orbits were averaged, and Pluto was ignored) showed that solar system itself was chaotic. Laskar suggested that secular resonances involving the terrestrial planets were responsible for the chaotic motion45 ,46 . Subsequent integrations of a complete model by Sussman and Wisdom15  confirmed that the full solar system was chaotic with $`T_L5`$ Myr. Sussman and Wisdom also found that the four giant planets by themselves appeared to form a chaotic system.

The chaos seen amongst the four giant planets arises from a three-body resonance involving Jupiter, Saturn, and Uranus16 . The orbital period of Uranus is nearly seven times that of Jupiter. In terms of the orbital frequency $`n=2\pi /T`$, the difference $`n_J7n_U`$ is equal to $`5n_S2n_J`$. The relevant resonant terms in the potential experienced by Uranus are proportional to $`M_JM_S/(25n_S/n_J)`$. Because $`(25n_S/n_J)1/85`$ appears downstairs, it is known as a “small denominator”; it enhances the potential experienced by Uranus by a factor of $`85`$. It was of course exactly this small denominator that Laplace used to explain the origin of the great inequality, and which is implicated in the removal of objects from the 2:1 resonance in the asteroid belt.

The width of a single component of this three-body resonance is tiny, comparable to the radius of Uranus, but so is the separation between components; the resonances just overlap. The predicted Lyapunov time is $`T_L10`$ million years. Because Uranus is so much less massive than either Jupiter or Saturn, its orbital angular momentum is substantially less than that of the two gas giants. Hence it is subject to the gambler’s ruin. The diffusion time (roughly the time before Uranus is ejected) is $`10^{18}`$ years, much longer than the current age of the universe. The theory also predicts the location and Lyapunov times of other chaotic zones near the present orbit of Uranus. Detailed numerical integrations verify the presence and Lyapunov times of all these zones, as well as illustrating the transitions between libration and rotation in the relevant resonances. Current computational resources are inadequate to test the predicted diffusion time.

Chaos among the terrestrial planets

As noted above, numerical integrations of the full solar system, inner planets included, show evidence of chaos with a Lyapunov time of $`T_L5`$ million years15 ,17 . Unlike the outer planets, the source of this chaos has not been convincingly established. Laskar pointed to what are called a “secular resonances” between the terrestrial planets as a candidate source of the chaos45 . He found what appears to be an alternation between circulation and libration in the angles $`\sigma _1(\varpi _1\varpi _5)(\mathrm{\Omega }_1\mathrm{\Omega }_2)`$ and $`\sigma _22(\varpi _4\varpi _3)(\mathrm{\Omega }_4\mathrm{\Omega }_3)`$. In these relations $`\varpi _1`$ refers to the orientation of Mercury’s apsidal line, which is the line from the sun to the point of Mercury’s orbit closet to the sun, while $`\mathrm{\Omega }_1`$ refers to the orientation of the Mercury’s nodal line; the latter is defined by the intersection of the orbital plane of Mercury with the orbital plane of the Earth. Similar definitions apply to the elements for Venus (2), Earth (3), Mars (4), and Jupiter (5). However, these two resonances do not interact directly with each other, so by themselves are they unlikely to produce large scale chaos. Later integrations46  identified a third resonance, $`\sigma _3(\varpi _4\varpi _3)(\mathrm{\Omega }_4\mathrm{\Omega }_3)`$. Laskar found that it librated while $`\sigma _2`$ rotated, and vice-versa. Since $`\sigma _2`$ and $`\sigma _3`$ involve the same degrees of freedom, they are a more promising candidate for overlapping resonances. Sussman and Wisdom, employing a more realistic (unaveraged) model, confirmed the alternate libration and rotation of $`\sigma _2`$, but never saw $`\sigma _3`$ librate; this does not rule out the possibility that the over lap of the associated resonances causes the chaos. However, their assessment that “no dynamical mechanism for the observed chaotic behavior of the solar system has been clearly demonstrated” seems warranted, at least for the terrestrial planets. Without a clear identification of the source of the chaos it is not possible to use an analytic development, such as was used for the outer planets, to confirm the Lyapunov time and then estimate the time scale for diffusion of the system.

Chaos in other planetary systems

By analogy with our solar system, multiplanet systems seem likely. In fact three Jupiter-mass objects orbit Upsilon Andromeda47 . Any multiple planet system is subject to the instabilities described in this paper. All of the known planetary systems orbiting solar type stars have ages of $`10^9`$ years or greater. This fact, together with the Copernican assumption that we are not observing the system at a privileged time, such as immediately after a recent planetary ejection, can be used to put limits on the mass and orbital elements of the planets, quantities which cannot be tightly constrained or even directly observed using current techniques48 ,49 . Dynamical constraints can also be placed on the existence of smaller (Earth mass) bodies in orbits near those of the Jupiter-mass objects. Observational searches for such low mass companions will have to await improved techniques.

Acknowledgements We thank Peter Goldreich for helpful conversations. We also thank the referees for their assistance in improving the manuscript. This research was supported by NSERC of Canada and NASA.

TEXT OF BOX 1 BEGINS HERE

Here we describe a simple model that captures many of the generic features of conservative or Hamiltonian dynamical systems. As stated earlier, at least two oscillators must be interacting or coupled for chaos to occur. We can add a second oscillator to the rigid pendulum by vertically (or horizontally) moving the pivot of the pendulum periodically. Much insight can be gained by carefully examining this and similarly simple models6 ,4 ,27 ,50

To describe the motion we now need two angles, say the angle of the pendulum, $`\theta `$, and the phase of the pivot, $`\psi `$. In addition we need two angle derivatives, $`\dot{\theta }`$ or $`p=\dot{\theta }/\beta `$, and $`\dot{\psi }`$. The Hamiltonian is

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2}}\beta p^2+\alpha I+\left[ϵ+2K\mathrm{cos}\psi \right]\mathrm{cos}\theta `$
$`=`$ $`{\displaystyle \frac{1}{2}}\beta p^2+\alpha I+ϵ\mathrm{cos}\theta +K\mathrm{cos}(\theta \psi )+K\mathrm{cos}(\theta +\psi ).`$
The first line presents an amplitude-modulated pendulum; the second presents three phase-modulated pendula. See ref. 27  for a detailed development. Several constants depend on the physical characteristics of the driven pendulum: $`\beta `$ is related to moment of inertia of the pendulum; $`K`$ is a proportional to the amplitude of the driving; and $`ϵ`$ is proportional to the gravitational acceleration. The constant $`\alpha `$ is the driving frequency or $`\dot{\psi }`$.

There are three resonances (cosine terms) which dominate the motion when their arguments are slowly varying: $`\dot{\theta }0`$, $`\dot{\theta }\dot{\psi }`$, and $`\dot{\theta }\dot{\psi }`$. We can visualize the trajectories by plotting $`p`$ versus $`\theta `$ whenever the driving reaches a specific phase, $`\alpha =0`$ for example. The resulting diagram is called a surface of section.

Figures 7a-7c show surfaces of section for $`\beta =1`$, $`ϵ=1`$, and $`K=0.3`$. The panels correspond to decreasing driving frequencies of $`\alpha =7`$, $`3`$, and $`0.3`$, respectively.

Here, the natural oscillation frequency of the pendulum is $`\sqrt{|ϵ|\beta }=1`$. The resonant islands are located at $`p=0`$ and $`\beta p=\pm \dot{\psi }=\pm \alpha `$. The center resonant island has a half-width of $`\mathrm{\Delta }P=2\sqrt{\frac{|ϵ|}{|\beta |}}=2`$ and the other two have half-widths of $`\mathrm{\Delta }P=2\sqrt{2K/\beta }1.5`$. So, for the highest driving frequency (top panel) the centers of the resonant islands are separated by more than their widths. Aside from the addition resonant islands, two new features appear. First, two chains of secondary islands can be seen. Second, the separatrix of the main island has broadened into a fuzzy zone. This zone was traced out by a single trajectory. Each of the separatrices shows such a chaotic zone, although they are too small to discern. Nevertheless, what appear to be invariant curves corresponding to regular trajectories can still be seen. That both regular and chaotic trajectories should appear on the same surface of section is a generic feature of a divided phase space18 ,19 ,27 .

In the middle panel the driving frequency has been reduced, $`\alpha =3`$, such that now the islands, ignoring their interaction, would just touch each other. As the driving frequency approaches the natural oscillation frequency of the pendulum, the small chaotic zone seen in the top panel envelopes the three resonant islands. This is the result of the beginning of resonance overlap27 .

In the bottom panel the driving frequency has been further reduced, $`\alpha =0.3`$, such that now the resonances have widths much larger than their separations. The motion is chaotic near the separatrix, which appears as a single thick band surrounding a stable island, and is surrounded by invariant curves. As can be seen from the Hamiltonian, the three resonant terms can also be viewed as a single resonance with an oscillating width, $`\mathrm{\Delta }P=2\sqrt{\frac{|ϵ+2K\mathrm{cos}\psi |}{\beta }}`$. The width ranges from $`\mathrm{\Delta }P1.3`$ to $`\mathrm{\Delta }P2.5`$, the chosen values of $`ϵ`$ and $`K`$. Here, the chaotic zone is roughly that region over which the separatrix sweeps. This is termed modulational chaos51  or adiabatic chaos4 .

In discussing the origin of chaotic motion in the main text, we have treated the region of overlap as though it had no structure. This simplification is justified in some cases (for example in high order resonances in the asteroid belt) and for some applications (for calculating Lyapunov times), but is not always adequate.

For example, Wisdom noted that the eccentricity of asteroids in the 3:1 resonance could librate around two different centers, one at the classical “forced eccentricity” $`e_f0.05`$ and a second at $`e0.15`$. The chaotic region surrounded the separatrix between these two points, allowing asteroids to pass from low eccentricities ($`0.05`$) to high $`0.3`$, where they could be removed by encounters with Mars. Later it was discovered that there was a third libration center at even higher eccentricity, which was connected to the first two by chaotic orbits52 . This allowed asteroids to reach $`e>0.6`$; such objects are removed from the resonance, typically by plunging into the sun. This motion is made possible by the resonance overlap, but describing it naturally requires knowledge of the underlying phase space structure; alternately one can perform numerical simulations.

Direct brute force numerical integration indicates that plunging into the sun is the most likely fate of bodies injected into the 3:1 resonance, with the remainder being pushed beyond Saturn53 . These integrations find a removal time of about one million years. There is also direct observational evidence that this interpretation is correct because the boundaries of the chaotic zone closely match the boundaries to observed distribution of asteroids3 . Similar statements can be made about the other Kirkwood gaps and gaps in the outer asteroid belt6  (see figure 6).

A more realistic model for the 3:1 mean motion resonance, still for a single perturbing planet, would include terms of order $`I^2`$, and allow for (different) $`I`$ dependences in the factors $`K`$ multiplying the cosine terms involving $`\psi `$. Such a model reproduces the two libration centers in the $`\psi ,I`$ motion found by Wisdom. Because the analytic model that we, and Wisdom, employ is of low order in $`e`$, it misses the third libration center found in reference 25.

Still more complete dynamical models, in particular those that include the effect of Saturn on Jupiter’s orbit, show that secular resonances, i.e., resonances between the precession periods of the asteroid and those of Jupiter and Saturn overlap inside low order mean motion resonance such as the 2:154 ,55 . This overlap can enhance the rate of diffusion. Even this is not the entire story; Jupiter and Saturn are near (but not in) a 5:2 resonance. This “great inequality” enhances the rate of diffusion of asteroids in the 2:1 mean motion resonance with Jupiter56 .

```
Figure 7a caption:

Surface of section for the driven pendulum with rapid driving much
faster than the natural pendulum oscillation frequency.
Top panel of figure in Box 1.

Figure 7b caption:

Surface of section for the driven pendulum with driving frequency
approaching the natural pendulum oscillation frequency.
Middle panel of figure in Box 1.

Figure 7c caption:
Surface of section for the driven pendulum with driving frequency
less than the natural pendulum oscillation frequency.
Bottom panel of figure in Box 1.

```

TEXT OF BOX 1 ENDS HERE