text:
# Electron Energy Spread in Laser Cooling Process

## I Introduction

The interaction of a high energy electron beam with an intense laser flash is now considered in such projects as compact x-ray sources, $`\gamma \gamma `$ colliders, diagnostics of low emittance electron beams, laser cooling and others Angelo . The laser flash intensity in these projects may be so high that an electron can undergo several successive collisions passing through a photon bunch Telnov-83 ; Telnov-95 . It means the number of scattered photons is larger than one. Taking into account the discrete character of the Compton back-scattering process, one can see an analogy between the electron passage through a laser flash and through a condensed matter. The stochastic theory of particles penetration through a matter can be found in the book We . The model of multiple bremsstrahlung process which leads to emission of a few photons by a single electron was developed in brems . Here we consider the multiple Compton back-scattering process using the same approach. For simplicity we use the approximation of an uniform photon concentration in the laser flash:

$$n_L=\frac{A}{\omega _0}\frac{1}{\pi r_L^2l_L}=\mathrm{const}$$

where $`A`$ is the total energy of the laser flash; $`\omega _0`$ is the photon energy, $`r_L`$ and $`l_L`$ are radius and length of the photon bunch. In other words our approximation corresponds to the short laser pulse ($`l_Lz_R`$, $`z_R`$ being the Rayleigh length).

The kinetic equations have been obtained for monoenergetic electrons and for electrons with incident energy spread. These equations are transformed into the equations for the distributions moments. The method for an approximate solution of the equations has been derived for the mean energy and for the variance of the energy distributions. It has been shown that the energy distribution variance as a function of the target thickness is a curve with a maximum if the incident energy spread is small.

It is shown that the equation for calculation of the variance can be used to study the variation of the energy spread in the laser-electron storage ring (LESR) where an electron beam repeatedly encounters the laser pulses. It has been shown that the variation of the energy spread in this case depends on the incident value of the variance. The energy spread damping takes place if the incident energy spread is larger than some limit value.

The results of analytical calculations are compared with the data of Monte Carlo simulation.

## II Transport equations

Let us consider electrons with an incident energy $`\epsilon _0mc^2`$ travelling through a bunch of photons with an energy $`\omega _0`$. The process of an electron transport in a light target is a random process where both the collisions number along the path $`l`$ and the energy loss in each collision are random quantities. Therefore, the energy of electron $`\epsilon `$ passing a laser flash is random quantity too.

The angular deflection of an electron in the Compton back-scattering process is small ($`{\displaystyle \frac{2\omega _0}{mc^2}}`$) and may be neglected (the straight-ahead approximation). In this approximation the probability density function $`P_m(\epsilon |\epsilon _0,l)`$ (the subscript ”m” means ”monoenergetic”) describing the energy distribution of electrons travelling a path $`l`$ in the photon bunch obeys the adjoint balance equation (the Kolmogorov-Chapman equation) Feller ; We :

$`P_m(\epsilon |\epsilon _0,l)=(1s\mathrm{\Sigma }(\epsilon _0))P_m(\epsilon |\epsilon _0,ls)`$
$`+s\mathrm{\Sigma }(\epsilon _0){\displaystyle \underset{0}{\overset{\omega _{max}}{}}}{\displaystyle \frac{\mathrm{\Sigma }(\omega ;\epsilon _0)}{\mathrm{\Sigma }(\epsilon _0)}}P_m(\epsilon |\epsilon _0\omega ,l)𝑑\omega ,`$ (1)
where $`\mathrm{\Sigma }(\epsilon _0)`$ and $`\mathrm{\Sigma }(\omega ;\epsilon _0)`$ are the total and differential macroscopic cross-sections of the Compton scattering, $`s`$ is a small part of $`l`$, and $`\omega `$ is the energy of scattered photon $`(0\omega \omega _{max}),`$

$$\omega _{max}=\epsilon _0\frac{x}{1+x}$$

is the maximum value of the scattered photon energy,

$$x=\frac{4\omega _0\epsilon _0}{(mc^2)^2},$$

$$\mathrm{\Sigma }(\epsilon _0)=2n_L\sigma (\epsilon _0),$$

$$\mathrm{\Sigma }(\omega ;\epsilon _0)=2n_L\frac{d\sigma (\omega ;\epsilon _0)}{d\omega },$$

$`\sigma `$ and $`{\displaystyle \frac{d\sigma }{d\omega }}`$ are the total and differential cross-sections. $`n_L`$ is the concentration of laser photons in the bunch.

Note that $`\mathrm{\Sigma }(\epsilon _0)`$ is the mean number of collisions of an electron per unit path length and $`\mathrm{\Sigma }(\omega ;\epsilon _0)`$ is the mean number of an electron collisions with the energy loss in unit interval about $`\omega `$ per unit path length.

The first term in the right side of Eq. (1) corresponds to the electrons, which passes the path $`s`$ without collisions, and $`1s\mathrm{\Sigma }(\epsilon _0)`$ is corresponding probability. These electrons have to lose the energy $`\epsilon _0\epsilon `$ along the rest path $`ls`$.

The second term corresponds to the electrons, which undergo the first scattering passing the path $`s`$, and $`s\mathrm{\Sigma }(\epsilon _0)`$ is the corresponding probability. The energy of the electron after the first scattering equals $`\epsilon _0\omega `$, where $`\omega `$ is a random energy of the scattered photon, and $`\mathrm{\Sigma }(\omega ;\epsilon _0)/\mathrm{\Sigma }(\epsilon _0)`$ is the probability density function of $`\omega `$. In the limit $`s0`$ Eq. (1) gives the adjoint integro-differential equation for the function $`P_m(\epsilon |\epsilon _0,l)`$ brems ; stochastics :

$`{\displaystyle \frac{}{l}}P_m(\epsilon |\epsilon _0,l)+\mathrm{\Sigma }(\epsilon _0)P_m(\epsilon |\epsilon _0,l)`$
$`{\displaystyle \underset{0}{\overset{\omega _{max}}{}}}\mathrm{\Sigma }(\omega ;\epsilon _0)P_m(\epsilon |\epsilon _0\omega ,l)𝑑\omega =0`$ (2)
with boundary condition

$$P_m(\epsilon |\epsilon _0,l)|_{l=0}=\delta (\epsilon \epsilon _0),$$

$`\delta (\epsilon \epsilon _0)`$ being the Dirac $`\delta `$-function.

If the electrons incident on the light target have an energy distribution $`P_0(\epsilon )`$, the energy distribution of electrons behind the bunch is the convolution of $`P_0(\epsilon )`$ and $`P_m(\epsilon |\epsilon _0,l)`$:

$$P_c(\epsilon |l)=\underset{\epsilon }{\overset{\mathrm{}}{}}P_0(\epsilon ^{})P_m(\epsilon |\epsilon ^{},l)𝑑\epsilon ^{}$$
(3)
(the subscript ”c” means ”convolution”).

Eq. (2) can be transformed into the equation for the distribution moments

$$\overline{\epsilon _m^k}(\epsilon _0,l)=\underset{mc^2}{\overset{\epsilon _0}{}}\epsilon ^kP_m(\epsilon |\epsilon _0,l)𝑑\epsilon .$$

The equation for $`\overline{\epsilon _m^k}(\epsilon _0,l)`$ stochastics  is

$`{\displaystyle \frac{}{l}}\overline{\epsilon _m^k}(\epsilon _0,l)+\mathrm{\Sigma }(\epsilon _0)\overline{\epsilon _m^k}(\epsilon _0,l)`$
$`{\displaystyle \underset{0}{\overset{\omega _{max}}{}}}\mathrm{\Sigma }(\omega ;\epsilon _0)\overline{\epsilon _m^k}(\epsilon _0\omega ,l)𝑑\omega =0.`$ (4)
The boundary condition for the moments is $`\overline{\epsilon _m^k}(\epsilon _0,l)`$ is

$$\overline{\epsilon _m^k}(\epsilon _0,l)|_{l=0}=\epsilon _0^k.$$

The equation for the moments of the distribution $`P_c(\epsilon |l)`$ describing transport of electrons with an incident energy spread can be obtained from Eq. (3):

$$\overline{\epsilon _c^k}(l)=\underset{mc^2}{\overset{\mathrm{}}{}}\epsilon ^kP_c(\epsilon |l)𝑑\epsilon =\underset{mc^2}{\overset{\mathrm{}}{}}P_0(\epsilon ^{})\overline{\epsilon _m^k}(\epsilon ^{},l)𝑑\epsilon ^{}.$$
(5)
It is seen that they are expressed in terms of the moments for the monoenergetic electrons and the incident energy distribution.

If the energy loss of an electron in one collision is small, the integro-differential equation for the moments (4) can be transformed by the Taylor expansion of integrands:

$`\overline{\epsilon _m^k}(\epsilon _0\omega ,l)`$ $``$ $`\overline{\epsilon _m^k}(\epsilon _0,l)\omega {\displaystyle \frac{}{\epsilon _0}}\overline{\epsilon _m^k}(\epsilon _0,l)`$
$`+{\displaystyle \frac{1}{2}}\omega ^2{\displaystyle \frac{^2}{\epsilon _0^2}}\overline{\epsilon _m^k}(\epsilon _0,l).`$
This gives the partial differential equation:

$`{\displaystyle \frac{}{l}}\overline{\epsilon _m^k}(\epsilon _0,l)`$ $`+`$ $`\beta (\epsilon _0){\displaystyle \frac{}{\epsilon _0}}\overline{\epsilon _m^k}(\epsilon _0,l)`$ (6)
$``$ $`{\displaystyle \frac{1}{2}}\gamma (\epsilon _0){\displaystyle \frac{^2}{\epsilon _0^2}}\overline{\epsilon _m^k}(\epsilon _0,l)=0.`$
The quantities $`\beta (\epsilon _0)`$, and $`\gamma (\epsilon _0)`$ in (6) are the moments of the macroscopic differential cross-section:

$$\beta (\epsilon _0)=\underset{0}{\overset{\omega _{max}}{}}\omega \mathrm{\Sigma }(\omega ;\epsilon _0)𝑑\omega ,$$

$$\gamma (\epsilon _0)=\underset{0}{\overset{\omega _{max}}{}}\omega ^2\mathrm{\Sigma }(\omega ;\epsilon _0)𝑑\omega .$$

Two first terms in Eq. (6) correspond to the continuous slowing down approximation, whereas the third term gives approximate description of the energy loss fluctuations.

The Eq. (6) written for $`k=1,2`$ can be transformed into the equation for the variance

$$\sigma _m^2(\epsilon _0,l)=\overline{\epsilon _m^2}(\epsilon _0,l)\overline{\epsilon }_m^2(\epsilon _0,l).$$

The equation has a form stochastics

$`{\displaystyle \frac{}{l}}\sigma _m^2(\epsilon _0,l)+\beta (\epsilon _0){\displaystyle \frac{}{\epsilon _0}}\sigma _m^2(\epsilon _0,l){\displaystyle \frac{1}{2}}\gamma (\epsilon _0){\displaystyle \frac{^2}{\epsilon _0^2}}\sigma _m^2(\epsilon _0,l)`$
$`=\gamma (\epsilon _0)\left({\displaystyle \frac{}{\epsilon _0}}\overline{\epsilon _m}(\epsilon _0,l)\right)^2.`$ (7)
The Eq. (5) for the energy distribution moments in the case of the electron beam with an incident energy spread can be simplified if the distribution $`P_0(\epsilon )`$ is a symmetric function with maximum at a point $`\epsilon _0`$ and small variance

$$\sigma _0^2=\underset{\mathrm{}}{\overset{\mathrm{}}{}}(\epsilon \epsilon _0)^2P_0(\epsilon )𝑑\epsilon .$$

In this case the Taylor expansion of the function $`\overline{\epsilon _m^k}(\epsilon ,l)`$ in (5) makes it possible to express the moments of the distribution $`P_c(\epsilon |l)`$ in terms of the moments $`\overline{\epsilon _m^k}(\epsilon _0,l)`$ corresponding to monoenergetic electrons:

$$\overline{\epsilon _c^k}(l)=\overline{\epsilon _m^k}(\epsilon _0,l)+\frac{\sigma _0^2}{2}\frac{^2}{\epsilon _0^2}\overline{\epsilon _m^k}(\epsilon _0,l).$$
(8)
The equations (8) written for $`k=1,2`$ can be easily transformed into the equation for the variance

$$\sigma _c^2(l)=\overline{\epsilon _c^2}(l)\overline{\epsilon _c}^2(l).$$

The equation is

$`\sigma _c^2(l)`$ $`=`$ $`\sigma _m^2(\epsilon _0,l)`$ (9)
$`+{\displaystyle \frac{\sigma _0^2}{2}}\left({\displaystyle \frac{^2}{\epsilon _0^2}}\sigma _m^2(\epsilon _0,l)+2({\displaystyle \frac{}{\epsilon _0}}\overline{\epsilon _m}(\epsilon _0,l))^2\right).`$
It describes the variance of the energy distribution for the electron beam with an incident variance $`\sigma _0`$ in terms of stochastic characteristics for monoenergetic electrons.

## III Compton scattering cross-sections and interaction coefficients

In this paper we restrict our consideration to the linear Compton scattering. The differential cross-section of this process for relativistic electron is

$`{\displaystyle \frac{d\sigma (y;x)}{dy}}={\displaystyle \frac{2\pi r_0^2}{x}}`$
$`\times \left(1y+{\displaystyle \frac{1}{1y}}{\displaystyle \frac{4y}{x(1y)}}+{\displaystyle \frac{4y^2}{x^2(1y)^2}}\right),`$ (10)
where $`y={\displaystyle \frac{\omega }{ϵ_0}}`$ and $`r_0={\displaystyle \frac{e^2}{mc^2}}`$ is the classical radius of electron Telnov-83 .

In the energy region of our interest the invariant dimensionless parameter $`x`$ is small $`(x1)`$ and the interaction coefficients $`\beta (\epsilon _0)`$, and $`\gamma (\epsilon _0)`$ are described by the approximate formulas

$$\beta (\epsilon _0)\frac{\mathrm{\Sigma }_0}{2}\epsilon _0x$$
(11)
$$\gamma (\epsilon _0)\frac{7}{20}\mathrm{\Sigma }_0\epsilon _0^2x^2,$$
(12)
where $`\mathrm{\Sigma }_0=2n_L\sigma _T`$, $`\sigma _T={\displaystyle \frac{8}{3}}\pi r_0^2`$ being the Thomson cross-section.

The quantities $`\beta (\epsilon _0)`$ and $`\gamma (\epsilon _0)`$ are the mean energy loss and the mean squared energy loss of an electron per unit path length.

## IV Solution of equations for monoenergetic electrons

The partial differential equations (6), (7) with the interaction coefficients (11), (12) describe the mean energy and the variance of the energy distribution for incidentally monoenergetic electrons passing a path $`l`$. The equations can be solved by the iteration method. In the first approximation, where the terms with the second derivatives are neglected, the solutions are

$$\overline{\epsilon _m}(\epsilon _0,l)=\frac{\epsilon _0}{1+\frac{1}{2}\overline{n}x}$$
(13)
$$\sigma _m^2(\epsilon _0,l)=\frac{7}{20}\frac{\epsilon _0^2x^2\overline{n}}{\left(1+\frac{1}{2}\overline{n}x\right)^4},$$
(14)
where $`\overline{n}=\mathrm{\Sigma }_0l`$ is the mean number of collisions in the light target stochastics .

The mean value (13) is a monotonically decreasing function of variable $`l`$. The Taylor expansions show that for small $`l`$ $`(\overline{n}x1)`$

$$\overline{\epsilon _m}(\epsilon _0,l)\epsilon _0(1\frac{\overline{n}x}{2})$$

and for large $`l`$ $`(\overline{n}x1)`$

$$\overline{\epsilon _m}(\epsilon _0,l)\frac{2\epsilon _0}{\overline{n}x}\left(1\frac{2}{\overline{n}x}\right).$$

Similarly, it can be shown that the variance (14) grows with $`l`$ for thin targets:

$$\sigma _m^2(\epsilon _0,l)\frac{7}{20}\overline{n}\epsilon _0^2x^2(12\overline{n}x)$$

and decreases for large $`l`$:

$$\sigma _m^2(\epsilon _0,l)\frac{28}{5}\frac{\epsilon _0^2}{\overline{n}^3x^2}(1\frac{8}{\overline{n}x}).$$

It follows from (14) that the variance $`\sigma _m^2(\epsilon _0,l)`$ has a maximum at the point, where

$$\overline{n}=\frac{2}{3}\frac{1}{x}$$

or

$$\frac{\overline{\epsilon }_m(\epsilon _0,l)}{\epsilon _0}=\frac{3}{4}.$$

The energy spread damping for large $`l`$ is due to the energy loss $`\beta (\epsilon )`$ per unit path length is proportional to $`\epsilon ^2`$ ,therefore, higher energy electrons in a bunch lose more energy than lower energy electrons Telnov-97 ; Huang1 .

Using Eqs. (14) and (13) one can derive the formula

$$\sigma _m^2(\epsilon _0,l)=\frac{7}{10}x\frac{\overline{\epsilon }}{\epsilon _0}\left(1\frac{\overline{\epsilon }}{\epsilon _0}\right),$$

which was earlier obtained in Telnov-97 .

## V Solution of equations for electron beam with incident energy spread

Substitution of Eqs. (13), (14) in Eqs. (8) and (9) makes it possible to calculate the mean energy and the variance of the energy distribution for the beam with the incident energy spread $`\sigma _0`$.

These calculations show that the mean value $`\overline{\epsilon _c}(l)`$ is decreased monotonically with $`l`$:

$$\overline{\epsilon _c}(l)\epsilon _0\left(1\frac{\overline{n}x}{2}(1+\frac{\sigma _0^2}{\epsilon _0^2})\right)$$

for small $`l`$ and

$$\overline{\epsilon _c}(l)\frac{2\epsilon _0}{\overline{n}x}\left(1\frac{2}{\overline{n}x}(1+\frac{\sigma _0^2}{\epsilon _0^2})\right)$$

for large $`l`$.

But the dependence of variance $`\sigma _c^2(l)`$ on $`l`$ is determined by the incident energy spread. The Taylor expansion of Eq. (9) in powers of $`l`$ gives

$$\sigma _c^2(\epsilon _0,l)\sigma _0^2+\frac{7}{20}\epsilon _0^2\overline{n}x^2\left(1\frac{40\sigma _0^2}{7\epsilon _0^2x}(1\frac{21x}{20})\right),$$

whereas for large $`l`$ it is decreased with $`l`$:

$$\sigma _c^2(\epsilon _0,l)\frac{28}{5}\frac{\epsilon _0^2}{\overline{n}^3x^2}\left(1\frac{8}{\overline{n}x}\left(1+\frac{\sigma _0^2}{\epsilon _0^2}(1\frac{10}{7x})\right)\right).$$

It means that the variance $`\sigma _c^2(l)`$ has a maximum at some $`l`$ if $`\sigma _0^2\sigma _{max}^2`$, where

$$\sigma _{max}^2=\frac{7}{40}x\epsilon _0^2\frac{1}{1\frac{21}{20}x},$$
(15)
For $`\sigma _0=\sigma _{max}`$ this maximum is at the point $`l=0`$. If $`\sigma _0>\sigma _{max}`$ the variance $`\sigma _c^2(l)`$ is a decreasing function of $`l`$.

## VI Multiple crossing of laser bunch

It was pointed out in Huang1 ; Huang2  that effective radiative cooling of an electron beam can be realized in LESR where electrons repeatedly interact with an intense laser pulses. Eq. (9) can be used to study the variation of the energy spread for electrons multiply crossing the light target. If the lost energy is restored by an rf accelerating system after each turn, the variance of energy distribution after n-th turn can be calculated using Eq. (9) with $`\sigma _0^2`$ replaced by the variance after (n-1)-th turn:

$$\sigma ^2(n)=\sigma _m^2(\epsilon _0,l)+\eta \sigma ^2(n1),$$
(16)
where

$$\eta =\frac{1}{2}\left(\frac{^2}{\epsilon _0^2}\sigma _m^2(\epsilon _0,l)+2(\frac{}{\epsilon _0}\overline{\epsilon _m}(\epsilon _0,l))^2\right).$$
(17)
It follows from (16) that for $`|\eta |1`$

$$\sigma ^2(\mathrm{})=\sigma _m^2(\epsilon _0,l)(1+\eta +\eta ^2+\eta ^3+\mathrm{})=\frac{\sigma _m^2(\epsilon _0,l)}{1\eta }.$$
(18)
Let us point out that the multiple crossing of the laser bunch leads to the energy spread damping if the incident energy spread $`\sigma _0^2>\sigma ^2(\mathrm{})`$. In the opposite case ($`\sigma _0^2<\sigma ^2(\mathrm{})`$) the energy spread is increased.

Limit value of energy spread $`\sigma ^2(\mathrm{})`$ depends on the thickness of laser bunch especially for high energy electrons. For small $`l`$

$$\sigma ^2(\mathrm{})\sigma _{max}^2\left(1\overline{n}x\frac{(15+28x)}{(2021x)}\right),$$

whereas for large $`l`$

$$\sigma ^2(\mathrm{})\frac{28}{5}\frac{\epsilon _0^2}{\overline{n}^3x^2}.$$

It is easy to see that for small $`x`$

$$\frac{\sigma (\mathrm{})}{\epsilon _0}\frac{\sigma _{max}^2}{\epsilon _0}\sqrt{\frac{7}{40}x}=\sqrt{\frac{7}{10}\gamma \frac{\lambda _e}{\lambda _L}},$$

where $`\lambda _L`$ is the laser radiation wavelength and $`\lambda _e={\displaystyle \frac{h}{mc}}`$ is the Compton wavelength of electron. This result coincides with that one given in Huang2 . But our approach allows to estimate the number of turns in a LESR to reach the limit value of the energy spread.

## VII Numerical results

Numerical calculation of the mean energy and the variance of the energy distribution were made for $`100MeV`$ electrons interacting with $`1.24eV`$ laser photons and for $`5GeV`$ electrons interacting with $`2.48eV`$ photons. Eqs. (8) and (9) were used to calculate the relative energy spread $`{\displaystyle \frac{\sigma _c(l)}{\overline{\epsilon _c}(l)}}`$. The results are given in Fig. 1 for 3 values of the incident energy spread. They agree with the data obtained using Eq. (15) of Telnov3 .

It is seen from the figure that the relative energy spread of the electron beam as a function of the collisions number $`\overline{n}`$ is a decreasing function if the incident energy spread is large enough. But for the beams with small incident spread the ”heating” of the beam takes place instead of ”cooling” for small $`l`$.

The dependence of the energy distribution variance on the number of turns in LESR is shown in Fig. 2 for two incident energy spreads. The calculations were made for LESR suggested for laser beam cooling in Huang2 . In this case $`\omega _0=1.24eV,\epsilon _0=100MeV`$, the mean energy of scattered photon equals $`\overline{\omega }={\displaystyle \frac{1}{2}}\epsilon _0x=95keV`$, the average ring radius equals $`1m`$, and the energy loss per turn $`\mathrm{\Delta }\omega =25keV`$. The mean number of collisions in laser bunch $`\overline{n}=\frac{\mathrm{\Delta }\omega }{\overline{\omega }}=0.27`$. It is seen from the figure that the asymptotic value of the energy spread is established in $`510^3`$ turns and it takes $`100\mu sec.`$

It is seen from Fig. 2 that the energy spread of electrons multiply crossing the light target is decreased or increased depending on its incident value.

The limit value of the relative energy spread $`\sigma (\mathrm{})/\epsilon _0`$ as a function of the light target thickness is given in Fig. 3. It is seen from the figure that the variance of the energy distribution is decreased with the target thickness increasing.

## VIII Monte Carlo simulation

It is supposed in our Monte Carlo simulation of the back Compton scattering that the incident energy distribution of electrons is Gaussian with given variance $`\sigma _0^2`$. The number of electron collisions with laser photons is supposed to be random. It is selected from the Poisson distribution with fixed $`\overline{n}`$. The simulation of individual collisions is carried out in the electron rest frame using the Klein and Nishina formula with the Lorentz transformation to the lab system. In the same way as in analytical calculations above we neglected the angular deflection of electrons but accounted for the energy decreasing after each collision.

Fig. 4 shows the energy distributions of $`5GeV`$ electrons for several values of the light target thickness. It should be pointed out the discontinuous of the spectra for small $`\overline{n}`$ due to the single scattered electrons. It is seen from the figure that the energy spread increases with the target thickness increasing for small $`\overline{n}`$ and decreases for such $`\overline{n}`$, which are greater than $`\sigma _{max}^2`$ determined by Eq. (15).

Fig. 5 and 6 shows the electron energy distributions in LESR for several number of turns.

All results are obtained with statistics more than $`10^6`$ trajectories.

It should be pointed out that the analytical formulas obtained in the continuous slowing down approximation with approximate consideration of the energy loss fluctuations are inapplicable for highly relativistic electrons if the mean energy loss in one collision is comparable to $`\epsilon _0`$. The Monte Carlo method has to be used in this case.

Fig. 7 shows the relative energy spread as a function of the light target thickness calculated in continuous slowing down approximation and by the Monte Carlo method for $`250GeV`$ electrons. Fig. 8 shows the energy spectra of these electrons. Notice that the Monte Carlo technique allows to get the results with taking into account the energy dependence of the interaction cross-sections and the lateral distribution of photons in the laser flash.


format: text

text:
# The Late-Time Radio spectrum of SN1993J

## 1 Introduction

The most widely accepted model for the radio emission from supernovae is the standard interaction model, in which the observed radio (synchrotron) emission is associated with the interaction region between the supernova ejecta and the pre-supernova wind of the progenitor star. The radio emission shows an early rise and a late decline due to a low-frequency absorption process. Chevalier (che82 (1982)) considered a number of mechanisms that could account for the absorption: (i) free-free absorption by material in front of the forward shock front, (ii) free-free absorption by material within the emitting region, (iii) synchrotron self-absorption, and (iv) the Razin-Tsytovich effect. At the beginning of the 1980’s, the available radio light curves could be well explained by synchrotron radio emission that was partially free-free absorbed by ionized thermal material in the circumstellar wind of the progenitor star (e.g. Weiler et al. wei86 (1986)). The availability of more complete light curves for these supernovae (e.g. Weiler et al. wei89 (1989), wei90 (1990)) provided further support for this mechanism. However, as new radio supernovae have been detected, the situation has become more complex. For example, the radio light curves of SN1986J showed a different shape in the rise to maximum, which was interpreted by Weiler et al. (wei90 (1990)) as evidence for internal free-free absorption. These authors suggested that a filamentary, clumpy-like structure of thermal absorbers could also be the source of the slow radio turn-on. Chevalier (che98 (1998)) has recently argued that synchrotron self-absorption is likely to be the dominant absorption mechanism for some radio supernovae, in particular for SN1987A, and for Type Ib and Type Ic supernovae. Finally, observations of SN1993J seem to indicate that more than one absorption mechanism plays a role (Fransson and Björnsson fb98 (1998); Pérez-Torres et al. pam01 (2001); Mioduszewski et al. mio01 (2001)).

SN1993J is to date the radio supernova whose evolution has been monitored in greatest detail. The brightness and proximity of SN1993J in M81 offered an unprecedented occasion for VLBI studies, almost since the supernova explosion on 28 March 1993. For the first time, a shell-like structure in a young radio supernova was discovered (Marcaide et al. 1995a ), in agreement with the circumstellar interaction model. Marcaide et al. (1995b ) showed that the expansion of SN1993J was self-similar, and produced the first movie of an expanding supernova. Later on, Marcaide et al. (mar97 (1997)) reported on the deceleration of the SN1993J expansion, a result recently confirmed by Bartel et al. (bartel00 (2000)). The modeling of the radio light curves of SN1993J has also given deeper insight into the supernova phenomenon. It seems now clear that the progenitor of SN1993J had a substantial mass-loss rate during the late stages of its evolution ($`\dot{M}`$$`10^5`$ $`M_{}\mathrm{yr}^1`$, e.g. Van Dyk et al. vdyk94 (1994)). The power-law density profile of the circumstellar medium around SN1993J seems to be different from the standard one (Lundqvist lun94 (1994); Van Dyk et al. vdyk94 (1994); Marcaide et al. mar97 (1997); Immler et al. immler01 (2001); Mioduszewski et al. mio01 (2001)), although Fransson and Björnsson (fb98 (1998)) have succeeded in modeling the SN1993J radio light curves with an $`s`$=2 density profile for the CSM.

Based on numerical modelings of the available radio data for SN1993J, Fransson and Björnsson (fb98 (1998)) and Pérez-Torres et al. (pam01 (2001)) made predictions about its late-time, low frequency ($`<1.4`$ GHz) radio light behaviour. However, observational results of the low-frequency radio emission from SN1993J were lacking. Thus, we carried out VLA observations between $``$322 and $``$15000 MHz aimed at checking the above modeling efforts, as well as a way to discern the main absorbing processes acting in SN1993J.

## 2 Observations and Data Reduction

We observed SN1993J using 25 antennas of the VLA in its most extended configuration (A, with baseline lengths from 0.68 up to 36.4 km) at L (18 cm) and P (90 cm) bands from UT 13:55 to UT 15:25 on 17 December 2000, and at U (14.9 GHz), X(8.5 GHz), C(4.9 GHz), L(1.3-1.7 GHz), and P(0.32 GHz) bands from UT 08:40 to UT 11:40 on 21 December 2000. Each frequency band was split into two intermediate frequencies (IFs). Data processing was made using the Astronomical Image Processing System (AIPS). We first describe briefly the observations at the higher frequencies, and then emphasize and describe in detail those at lower frequencies, particularly those in the P-band, deemed to be a discriminator of the dominant absorption mechanism in SN1993J. Table 1 summarizes our observations. The flux densities reported were obtained by combining the data from both intermediate frequencies (IFs) at each frequency band, except at L-band, for which the frequency separation was large enough that we could use each flux density value independently.

### 2.1 High-frequency (U-, X-, and C-band) observations

We used a standard continuum mode at U, X, and C bands, thus covering a bandwidth of 50 MHz per IF and frequency band. At all three bands, we observed SN1993J phase-referenced to 0917+624. The source 0917+624 also served as the phase reference for the system calibration. At U- and X-band, 3C286 was used as primary flux calibrator (assumed of constant flux density), while at C-band 3C147 was used instead. At U-band, the integration times on SN1993J and 0917+624 were $`\tau _{93\mathrm{J}}`$=21 min and $`\tau _{0917}`$ = 7 min, while at X-band these were 11.5 and 6 min, and at C-band 12 and 7 min, respectively.

### 2.2 Low-frequency (L- and P-band) observations

Observations were done in spectral-line mode in order to detect the likely presence of radio frequency interference (RFI) in the data, and also to avoid bandwidth smearing. In addition, this observing mode allows a more careful editing. 3C286 was used at both bands as primary flux calibrator, as well as bandpass calibrator.

L-band observations. Each of the two 12.5-MHz IFs at L-band (centered at 1.34 and 1.67 MHz, respectively) were split into 8 spectral channels. Observations of SN1993J ($`\tau _{93\mathrm{J}}`$ = 11.5 min) were phase-referenced to 0917+624 ($`\tau _{0917}`$ = 10.5 min). 0917+624 was also used as secondary flux calibrator at L-band. Our final hybrid map at L-band (not shown here) was made at an effective observing frequency of 1.50 GHz. The peak of brightness of the map (97.4 mJy beam-1) corresponds to the nucleus of M81, and the flux density of SN1993J is 31.9 mJy. The rms of the image background is 0.3 mJy/beam. The flux density values reported at L-band in Table 1 correspond to each IF, 1.34 and 1.67 GHz, respectively. Since the center frequencies in L-band are quite separated, they can be used as independent values for the purpose of fitting a spectral index (see below). We note that, though the distance between SN1993J and M81 (2$`\mathrm{}\mathrm{\hspace{0.17em}48}\mathrm{}`$) is much larger than the dimensions of the restoring beam ($`1.84\times 0.99`$ arcsec), wide field mapping is necessary if one aims at accurate determinations of the flux densities of SN1993J at low frequencies.

P-band observations. At P-band, each of the two 3.125-MHz IFs (centered at 322 and 327 MHz, respectively) was split into 32 spectral channels, which permitted a detailed inspection of the data, looking for the presence of RFI. Since 0917+624 is a relatively weak source at P-band, unlike at higher frequencies, confusion of many strong sources in the field made it unsuitable as secondary calibrator. Therefore, the much stronger source J1206+6414 (3C268.3) was used instead. The total integration times were $`\tau _{93\mathrm{J}}`$78 min, and $`\tau _{3\mathrm{C}268.3}`$39 min. We imaged 3C268.3 to obtain an accurate value of its flux density. We noticed that the coordinates used for 3C268.3 at the correlator ($`\alpha =12^\mathrm{h}\mathrm{\hspace{0.17em}06}^\mathrm{m}\mathrm{\hspace{0.17em}24}\stackrel{s}{.}711`$, $`\delta =+64\mathrm{°}\mathrm{\hspace{0.17em}14}\mathrm{}\mathrm{\hspace{0.17em}40}\stackrel{}{.}939`$; J2000.0) were off from the actual ones by more than one arcmin. Indeed, the coordinates obtained using the AIPS task JMFIT were $`\alpha =12^\mathrm{h}\mathrm{\hspace{0.17em}06}^\mathrm{m}\mathrm{\hspace{0.17em}24}\stackrel{s}{.}8978`$, $`\delta =+64\mathrm{°}\mathrm{\hspace{0.17em}13}\mathrm{}\mathrm{\hspace{0.17em}36}\stackrel{}{.}425`$ (J2000.0). Such position offset for 3C268.3 translated into a similar one for SN 1993J and the whole sky around it. We corrected the phases of 3C268.3 for this offset.

With the 3C268.3 data corrected, we imaged SN1993J at P-band, aimed to determine its flux density. Since the primary beam of the VLA at P-band is so large ($`2\mathrm{°}\mathrm{\hspace{0.17em}20}\mathrm{}`$), bright sources far from the supernova could greatly affect our image. Further, to obtain a good signal-to-noise ratio (SNR) in the image, one is forced to map even distant, strong sources. The situation is similar to that encountered at L-band, where self-calibration was mandatory to improve the SNR of the image, but useless unless M81 was also included in the imaging process. At P-band, not only the M81 nucleus is important, but also far away, strong sources such as M82 need also to be taken into account to be able to apply self-calibration to the data. We used the AIPS task SETFC to create a list of fields with suspected sources in them, rather than making an enormous, computationally unmanageable image. It is worth noticing that the interferometric phases are initially so poorly calibrated that the reconstructed image will show very few sources, making very difficult the convergence of the hybrid mapping procedure. Therefore, we used the AIPS procedure FACES, which searchs for sources in the NVSS catalog, to supply an initial model of the sky, and thus improve the global phase calibration. We then applied wide-field imaging on this self-calibrated dataset: since the array configuration looks different to sources in different parts of the (large) primary beam, we handled this by computing the uv-coverage differently for many small fields within the task IMAGR. Sources above a threshold of 33 mJy (at a frequency of 1.4 GHz) were searched for as “a priori” candidates in the NVSS catalog, within a 2$`\times `$2 deg region around SN1993J. This resulted in 78 fields that were cleaned using wide-field mapping. Fig. 1 shows the inner 1$`\times `$1 deg of our final, composite image of the VLA-A field at P-band. For clarity, only contours above 20 mJy beam-1 have been drawn. Panel (b) shows the galaxy M 82, offset about $`0\stackrel{}{.}40`$ from SN1993J. The latter is shown in panel (c) along with the nucleus of its host galaxy M81. M82 has a flux of $``$10.5 Jy, at an epoch when SN1993J has a flux of a mere 71.1 mJy (see Table 1 for details of the observations).

## 3 The radio spectrum of SN1993J

We fitted the radio spectrum of SN1993J obtained on 17 and 21 December 2000 using three models for its radio emission. These models are described by the equations

$$S(\nu )=K_1\nu _5^\alpha $$
(1)
$$S(\nu )=K_1\nu _5^\alpha e^\tau $$
(2)
$$S(\nu )=K_1\nu _5^\alpha (1e^\tau ^{})/\tau ^{}$$
(3)
where $`\tau =K_2\nu _5^{2.1}`$, and $`\tau ^{}=K_3\nu _5^{2.1}`$. These equations represent a simple power-law spectrum (model 1), a power-law spectrum free-free absorbed by a screen of homogeneously distributed ionized gas (model 2), and a power-law spectrum free-free absorbed by a “clumpy” medium (model 3). To facilitate comparisons with the modeling of Van Dyk et al. (vdyk94 (1994)), we have used here their notation: $`\nu _5`$ is the frequency, in units of 5 GHz; $`\alpha `$ is the observed spectral index; $`K_1`$ is the flux density at 5 GHz, in mJy; and $`K_2`$ and $`K_3`$ are the optical depths at 5 GHz of the homogeneous and “clumpy” external absorbing media, respectively, at the epoch of our observations. Table 2 summarizes our results, and Fig. 2 shows the best fits to models 1 through 3. A simple power-law spectrum (model 1) yields a good fit to the data from 14.94 GHz down to 1.34 GHz, but departs significantly from the actual flux density at frequencies below 400 MHz, and can be ruled out. The best fits are given by models 2 (homogeneous absorbing medium) and 3 (clumpy absorbing medium), which yield an observed spectral index $`\alpha =0.67\pm 0.02`$. The frequency at which the free-free optical depth becomes unity is very similar for both model 2 ($`\nu _{\tau =1}=169`$ MHz) and model 3 ($`\nu _{\tau =1}=241`$ MHz). Note, however, that the very-low radio frequency behavior of the spectrum is significantly different (see Fig. 2). While model 2 (solid line) predicts a sharp fall-off below frequencies of $``$200 MHz, model 3 (dashed line) predicts a much milder flux density decrease. Unfortunately, the absence of radio measurements below 322 MHz does not allow us to rule out either of these models.

Thus, from Table 2 we cannot exclude, or prefer, any of the two alternative models.

To model the early-time ($`t<230`$ days) radio light curves of SN1993J, Van Dyk et al. (vdyk94 (1994)) invoked a combined “homogeneous-clumpy” model of the form $`S_\nu \nu _5^\alpha e^\tau (1e^\tau ^{})/\tau ^{}`$. For such model, we can also obtain acceptable solutions (e.g. $`\alpha =0.67;K_1=14.4;K_2=8\times \mathrm{\hspace{0.17em}10}^4;K_3=10^5`$) at the expense of increasing $`\chi _{\mathrm{red}}^2`$ by $`\mathrm{\hspace{0.17em}50}\%`$. If the fitting parameters are left unconstrained, this model yields an even better fit than models 2 and 3, but at the unaffordable price of a negative –hence unphysical– value for $`K_2`$, or $`K_3`$. Clearly, the combined ”homogeneous-clumpy” model cannot be preferred over models 2 or 3, but we also feel we cannot exclude it solely on the basis of a 50% increase of the reduced $`\chi ^2`$.

For a power-law relativistic electron distribution, $`N_E=N_0E^\gamma `$, the spectral index of the electron distribution $`\gamma `$ is related to the observed spectral index: $`\gamma =12\alpha `$. From our value of $`\alpha `$, it then follows that $`\gamma =2.34\pm 0.04`$. This value of $`\gamma `$ agrees well with that predicted by Fransson & Björnsson (fb98 (1998)) for $`t>1000`$ days. We should note that synchrotron self-absorption could play a role in the radio emission from SN1993J (Chevalier che98 (1998), Fransson & Björnsson fb98 (1998), Pérez-Torres et al. pam01 (2001), Mioduszewski et al. mio01 (2001)). In particular, the modeling of Fransson & Björnsson (see their Fig. 12), predicted a flux density of $`70`$ mJy at 0.32 GHz and for the approximate epoch of our observations, while the modeling of Pérez-Torres et al. predicted at that frequency a flux density larger than 100 mJy. Such an excess in flux of the latter modeling is likely related to the different spectral evolution obtained for the electron spectrum.

Fig. 3 shows the evolution of the observed spectral index of SN1993J from $``$70 up to 2820 days after the explosion. For each epoch, we fitted the available radio continuum data (Van Dyk, private communication) using model 2 as described above. While at early epochs the fits are not superb, they seem to be compatible with $`\alpha 1`$. However, for epochs $`>`$1000 days, there seems to be a clear trend: the observed spectral index, $`\alpha `$, becomes progressively less steep. This spectral evolution could be explained assuming significant energy losses, mainly associated with synchrotron, as well as Coulomb and expansion losses. Fransson and Björnsson (fb98 (1998)) proposed a model that reproduces reasonably well the observed spectral evolution of SN1993J. In their model, a constant fraction of the shocked thermal electrons –characterized by a constant spectral index– are injected and accelerated. These electrons lose their energy mainly due to synchrotron losses, thus steepening the integrated electron spectrum. Moreover, these authors predict that, for the physical parameters they are considering, radiative (synchrotron) losses start to become less important at $`t`$ 1000 days at frequencies below $``$2.3 GHz, whereas energy losses of the electrons due to adiabatic and Coulomb cooling contribute significantly below 2.3 GHz, which seems consistent with the observational results. We note that a multi-frequency monitoring of SN1993J at radio frequencies down to $`70`$ MHz, would likely prove to be very useful to further constrain the population of electrons responsible for the emission, as well as to better understand the late-time radio evolution of SN1993J.

Fig. 4 shows the fits of model 2 to all data points (solid line), and to all but the P-band data point (dashed line). Because any absorption model will obviously fit by turning down sharply below the lowest frequency data point, the availability of radio continuum measurements at the lowest possible frequencies is important, in particular as a radio supernova ages. Had the P-band data not been available, we could have wrongly concluded that the peak of the SN1993J spectrum is at $`\nu `$ 0.8 GHz at a level of $``$35 mJy. Having the P-data point available, we conclude that the peak is around, or below 0.3 GHz, and at least at a level of 70 mJy. (We note that a similar fit to all data, but using model 3 rather than model 2, would not change this result.) Similarly, the best fit to our SN1993J radio data using models 2 or 3 as described above, but not including the P-band measurement, yields $`\alpha =0.73`$, a somewhat steeper value than is obtained if all data points are used, $`\alpha =0.67`$, but within 2$`\sigma `$ of it.

## 4 Summary

We present VLA radio continuum measurements of SN1993J in M81 taken on 17 and 21 December, 2000, about 2820 days after the supernova explosion, spanning the frequency range from 0.32 to 14.9 GHz. We discuss in detail the P-band (322-327 MHz) observations, since the task of obtaining the flux density at these frequencies is much harder than at higher frequencies, due to confusion from many nearby and far away, strong sources, e.g. M82. We also point out the importance of having a radio spectral coverage that extends down to P-band as the supernova ages, since otherwise the conclusions can be wrong.

We discuss three models for the late radio emission of SN1993J: (i) a simple power-law spectrum, (ii) a power-law spectrum free-free absorbed by a screen of homogeneously distributed ionized gas, and (iii) a power-law spectrum free-free absorbed by a “clumpy” medium. We find that the best fit to the data is yielded by a power-law spectrum free-free absorbed by either a homogeneous, or a clumpy, distribution of ionized gas. However, a mixture of homogeneous and clumpy absorbers, while not preferred, cannot be totally ruled out. The radio spectrum between 0.32 and 14.94 GHz is well characterized by $`\alpha =0.67\pm 0.02`$ ($`S_\nu \nu ^\alpha `$), typical of an optically thin radio supernova, and is significantly different from that obtained at epochs $`<`$1000 days between 1.4 and 14.9 GHz, which indicates an ongoing spectral evolution in the radio emission from SN1993J. Since the spectral index $`\gamma `$ of the electron distribution ($`N_E=N_0E^\gamma `$), is related to the observed spectral index by $`\gamma =12\alpha `$, our value of $`\alpha `$ translates into $`\gamma =2.34\pm 0.04`$, which agrees well with that predicted by Fransson & Björnsson (fb98 (1998)) for $`t>1000`$ days. We also fit the available radio continuum data of SN1993J for the period from $``$70 up to 2820 days since the explosion, using a power-law spectrum free-free absorbed by a screen of homogeneously distributed ionized gas. The fit shows that the observed spectral index of SN1993J has been slowly evolving since $`t`$1000 days on, with $`\alpha `$ increasing from a value close to -1 to -0.67. This spectral evolution seems to suggest that, in addition to the radiative (synchrotron) losses, adiabatic cooling and ionization (Coulomb) losses at the lowest frequencies might be contributing significantly to the integrated electron spectrum. An extension of the low frequency observations further down to $``$70 MHz should help to fine tune the physical parameters of SN1993J, and discern between the homogeneous or clumpy nature of its circumstellar medium. From a technical point of view, such observations will be challenging.

###### Acknowledgements.

We thank Lucas Lara for his help with some aspects of the data reduction process, and Gianfranco Brunetti for discussions. We thank Kurt Weiler for making the VLA data publicly available. We are also grateful to an anonymous referee for a careful and constructive review of our manuscript. MAPT is grateful to the Instituto de Astrofísica de Andalucía for its hospitality during a short visit. This research has been supported by a Marie Curie Fellowship of the European Community (contract IHP-MCFI-99-1), and by the Spanish DGICYT grants AYA2001-2147-C02-01 and AYA2001-2147-C02-02. The VLA is an instrument of the National Radio Astronomy Observatory, which is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Incorporated.


text:
# A Collisional Family in the Classical Kuiper Belt

## 1 INTRODUCTION

According to linear secular theory, the inclination, $`i`$, and longitude of ascending node, $`\mathrm{\Omega }`$, of the orbit of a test particle embedded in a planetary system evolve with time, $`t`$, as

$`pi\mathrm{sin}\mathrm{\Omega }`$ $`=`$ $`p_{\mathrm{𝑓𝑜𝑟𝑐𝑒𝑑}}+p_{\mathrm{𝑓𝑟𝑒𝑒}}`$ (1)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}A_i\mathrm{sin}(f_it+\gamma _i)+i_{\mathrm{𝑓𝑟𝑒𝑒}}\mathrm{sin}(Bt+\gamma _{\mathrm{𝑓𝑟𝑒𝑒}})`$
$`qi\mathrm{cos}\mathrm{\Omega }`$ $`=`$ $`q_{\mathrm{𝑓𝑜𝑟𝑐𝑒𝑑}}+q_{\mathrm{𝑓𝑟𝑒𝑒}}`$ (2)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}A_i\mathrm{cos}(f_it+\gamma _i)+i_{\mathrm{𝑓𝑟𝑒𝑒}}\mathrm{cos}(Bt+\gamma _{\mathrm{𝑓𝑟𝑒𝑒}})`$
\[see, e.g., Chapter 7 of Murray & Dermott (1999)\]. Here $`f_i`$ and $`\gamma _i`$ are functions of the masses, orbital semi-major axes, inclinations, and nodes of the planets; $`A_i`$ depends on these same quantities and the semi-major axis of the particle; and $`B`$ is given by the planetary masses and orbital semi-major axes of the planets and of the particle. The number of eigenmodes describing the forced vertical motion equals $`N`$ and is usually equal to the number of planets in the system. The “free” or “proper” inclination, $`i_{\mathrm{𝑓𝑟𝑒𝑒}}`$, and $`\gamma _{\mathrm{𝑓𝑟𝑒𝑒}}`$ are constants specified by the initial $`i`$ and $`\mathrm{\Omega }`$. Analogous equations exist for the test particle’s eccentricity, $`e`$, and longitude of periastron, $`\stackrel{~}{\omega }`$.

Hirayama (1918) noted that certain asteroids share similar values of $`i_{\mathrm{𝑓𝑟𝑒𝑒}}`$, $`e_{\mathrm{𝑓𝑟𝑒𝑒}}`$, and semi-major axis, $`a`$. He termed these groups “families,” and christened the first three families Koronis, Eos, and Themis using the names of their respective first-discovered members. A family’s members are thought to be fragments of a collisionally disrupted parent body. Immediately following the disruption event, each fragment is characterized by osculating orbital elements $`a_0+\delta a`$, $`e_0+\delta e`$, and $`i_0+\delta i`$, where the subscript $`0`$ pertains to the parent body just prior to disruption. Fragments for which $`\delta a/a_0`$, $`\delta e/e_0`$, $`\delta i/i_01`$ have nearly the same $`i_{\mathrm{𝑓𝑟𝑒𝑒}}`$ and $`e_{\mathrm{𝑓𝑟𝑒𝑒}}`$. A non-zero $`\delta a`$ induces a change in the free precession frequency, $`\delta B`$. As the orbits of fragments differentially precess in their nodes and apses, their osculating $`i`$’s and $`e`$’s may span wide ranges, particularly if the magnitude of the forced term is close to that of the free term. However, the constants of integration, $`i_{\mathrm{𝑓𝑟𝑒𝑒}}`$ and $`e_{\mathrm{𝑓𝑟𝑒𝑒}}`$, will remain fixed barring non-secular effects (e.g., close encounters with planets, or mean-motion resonant effects), effectively recording the inceptive shattering event.

Classical Kuiper Belt Objects \[CKBOs; see, e.g., the reviews by Jewitt & Luu (2000), Farinella & Davis (2000), and Malhotra, Duncan, & Levison (2000)\], which we define to be those bodies having orbital semi-major axes between 43 and 47 AU, currently do not reside within low-order mean-motion resonances with the giant planets. Secular theory well describes the dynamical evolution of CKBOs on low-$`e`$, Neptune-avoiding orbits, modulo the effects of weak chaos (Torbett 1989; Torbett & Smoluchowski 1990; Duncan, Levison, & Budd 1995). Here we report on an elementary calculation of the free inclinations and free eccentricities of known CKBOs. We uncover evidence for the first collisional family in trans-Neptunian space. The candidate family is dynamically reminiscent of a sub-family of the Koronis population in the asteroid belt. Section 2 describes our computational method, and section 3 contains results and discussion.

## 2 METHOD

Osculating elements of CKBOs were downloaded from the Minor Planet Center (MPC; http://cfa-www.harvard.edu/iau/lists/TNOs.html). We selected only those objects observed at multiple apparitions and having fitted semi-major axes $`43<a(\mathrm{AU})<47`$. See Millis et al. (2002) for a discussion of how orbital uncertainties depend on the arclength of astrometric observations. The resultant list contained 125 CKBOs.

For each CKBO having an observed $`p_{\mathrm{𝑜𝑏𝑠}}i_{\mathrm{𝑜𝑏𝑠}}\mathrm{sin}\mathrm{\Omega }_{\mathrm{𝑜𝑏𝑠}}`$ and $`q_{\mathrm{𝑜𝑏𝑠}}i_{\mathrm{𝑜𝑏𝑠}}\mathrm{cos}\mathrm{\Omega }_{\mathrm{𝑜𝑏𝑠}}`$, we compute

$$i_{\mathrm{𝑓𝑟𝑒𝑒}}=\sqrt{(p_{\mathrm{𝑜𝑏𝑠}}p_{\mathrm{𝑓𝑜𝑟𝑐𝑒𝑑}})^2+(q_{\mathrm{𝑜𝑏𝑠}}q_{\mathrm{𝑓𝑜𝑟𝑐𝑒𝑑}})^2}$$
(3)
and similarly for $`e_{\mathrm{𝑓𝑟𝑒𝑒}}`$. To evaluate the forced terms, we employ the secular theory of Brouwer & van Woerkom (1950; hereafter BvK) for all eight planets of the solar system, as transcribed by Murray & Dermott (1999). Now the osculating elements of KBOs given by the MPC pertain to epochs between 1994 and 2002 A.D., and are referred to the ecliptic and equinox of J2000. Since the BvK solution defines $`t=0`$ to be 1900 A.D., we take $`t=100`$ yrs in equations (1) and (2) and ignore the small, 8-year spread in MPC epochs. Since the BvK solution refers to the ecliptic and equinox of B1950, we transform the MPC orbital angles ($`i`$, $`\mathrm{\Omega }`$, $`\stackrel{~}{\omega }`$) to the B1950 frame using formulae provided by Montenbruck (1989, pages 18–19). Differences between untransformed and transformed angles are no more than 1. We note also typographical errors in the transcription of the BvK solution by Murray & Dermott (1999); in their captions for Tables 7.2 and 7.3, $`e_{ij}`$ and $`I_{ij}`$ should instead read $`e_{ji}`$ and $`I_{ji}`$, respectively, in keeping with their notation throughout the rest of chapter 7 (subscript $`j`$ denotes the planet, while subscript $`i`$ denotes the eigenmode).

Our method for extracting proper elements is crude compared with modern techniques. We neglect terms higher than degree 2 in planetary eccentricities and inclinations; see Milani & Knežević (1994) to remedy this deficiency. We retain only terms that are linear in the eccentricity and inclination of the CKBO, forsaking the semi-numerical approach of Williams (1969) and more recent, purely numerical approaches (Knežević and Milani 2001) that avoid expansions of the test particle’s orbit altogether. Thus, we cannot expect accuracy in our results for highly inclined and eccentric CKBOs. Nonetheless, many observed CKBOs execute nearly circular and co-planar trajectories for which we expect the errors introduced by our analysis to be small. And for simplicity and ease of implementation, our procedure is probably bested only by Hirayama’s celebrated and successful (1918) analysis, which included only a truncated Jovian secular potential \[$`N`$ = 1 in equations (1) and (2)\]. It is in his pioneering spirit that we search for the first collisional family in the Kuiper Belt.

## 3 RESULTS AND DISCUSSION

Figure 1 displays the inclinations and eccentricities, observed and free, of CKBOs against their semi-major axes. While the observed elements betray no obvious clumping, this is not true for the free inclinations. We consider the population of KBOs having $`44a(\mathrm{AU})45`$ and $`1i_{\mathrm{𝑓𝑟𝑒𝑒}}(\mathrm{deg})3`$ to be a candidate collisional family. Table 1 provides the designations, semi-major axes, free elements, and $`H`$-magnitude inferred diameters of our 32 candidate family members.

The free inclinations of our candidate family cluster about $`2^{}`$, similar to the inclination of the invariable plane with respect to the ecliptic, $`i_{\mathrm{𝑖𝑛𝑣𝑎𝑟}}`$. This coincidence might lead one to suspect that these bodies represent primordial ones that coagulated in the invariable plane, rather than fragments generated by catastrophic disruption. In this alternative scenario, these bodies would appear in $`p`$-$`q`$ space at the point representing the inclination and node of the invariable plane referred to the ecliptic. Would all such bodies share the same $`i_{\mathrm{𝑓𝑟𝑒𝑒}}`$, without recourse to collisional genesis? There are two ways by which they might do so, but both paths are prohibitive. The first way is if $`i_{\mathrm{𝑓𝑜𝑟𝑐𝑒𝑑}}i_{\mathrm{𝑓𝑟𝑒𝑒}}`$ (see Figure 7.3 in Murray & Dermott), so that $`i_{\mathrm{𝑓𝑟𝑒𝑒}}i_{\mathrm{𝑖𝑛𝑣𝑎𝑟}}`$. The former condition fails to be satisfied; for every one of our candidate members, $`i_{\mathrm{𝑓𝑜𝑟𝑐𝑒𝑑}}i_{\mathrm{𝑓𝑟𝑒𝑒}}`$. The second option demands that all bodies appear in $`p`$-$`q`$ space over a timescale much shorter than the local forced precession timescale of $``$$`2\pi /f_82\times 10^5\mathrm{yr}`$. In this case, the forced inclination vectors of bodies appearing in $`p`$-$`q`$ space would be nearly identical since there would be insufficient time for the forced vectors to rotate. There would be differences in the amplitudes of the forced inclination vectors from body to body, but these would be small since $`\mathrm{\Delta }a/a`$ is small for our candidate family and the bodies are far from secular resonances. Consequently, the free inclination vectors of all bodies appearing at a single point in $`p`$-$`q`$ space would be nearly the same. However, for all members of our candidate family to have formed well within a time interval of $`2\times 10^5\mathrm{yr}`$ seems unreasonable. Current planetesimal coagulation calculations indicate KBO growth times of order $`10^7\mathrm{yr}`$ (Kenyon 2002, and references therein). By contrast, a collisional disruption event is virtually instantaneous; we proceed under the hypothesis of a collisional origin for these objects.

Our candidate family boundaries were chosen based on visual inspection of Figure 1c and the remarkable narrowing of the distribution of inclinations between panels 1a and 1c exhibited by our candidate family. A Kolmogorov-Smirnov (K-S) test comparing the distributions of observed and free inclinations for our candidate family yields a 0.1% probability that the two distributions are drawn from the same underlying distribution. We consider this probability sufficiently low that we believe the narrowing of the inclination distribution to be significant. No other subset of KBOs in Figure 1 evinces the same degree of sharpening when we transform from observed to free elements; we compute corresponding K-S probabilities of $``$20% for other subsets of KBOs.

Comparing the distribution of free inclinations of our candidate family to that of 22 KBOs having $`43a(\mathrm{AU})44`$ and $`1i_{\mathrm{𝑓𝑟𝑒𝑒}}(\mathrm{deg})4`$ reveals no significant difference; the K-S test yields a 38% probability that the two distributions are drawn from the same population. A comparison between our candidate family and the 12 KBOs having $`45a(\mathrm{AU})46`$ and $`1i_{\mathrm{𝑓𝑟𝑒𝑒}}(\mathrm{deg})3`$ gives a corresponding probability of 13%, again too high to claim a significant difference. Thus, our family, if real, might extend over the full gamut of semi-major axes from 43 to 46 AU. Subjectivity in the choice of family boundaries is a problem that afflicts even the relatively mature field of asteroidal dynamics, and clearly more objects would be desirable. We defer more sophisticated analyses of the significance of the clumping in $`i_{free}`$-$`a`$ space to future study. For the present, we reserve our attention to the 32 KBOs between $`a=44\mathrm{A}\mathrm{U}`$ and 45 AU that exhibit the sharpest concentration of free inclinations.

The free eccentricities of our candidate family are more broadly distributed than their free inclinations. This characteristic is reminiscent of the Koronis asteroid sub-family located at semi-major axes $`2.84a(\mathrm{AU})2.91`$. This sub-family exhibits an absolute width of $``$0.025 in $`e_{\mathrm{𝑓𝑟𝑒𝑒}}`$ and an absolute width of $``$0.005 in $`\mathrm{sin}i_{\mathrm{𝑓𝑟𝑒𝑒}}`$ (Bottke et al. 2001).111The Koronis family divides into two groups (“sub-families”), one located at $`2.84a(\mathrm{AU})2.91`$ and a second located at $`2.92a(\mathrm{AU})2.95`$. The free eccentricities in the latter sub-family are markedly higher than those in the former, reflecting excitation via passage through a secular resonance located at 2.92 AU (Bottke et al. 2001). The free eccentricities and free inclinations of the former sub-family are thought to be of collisional origin, and we employ this sub-family as the analogue of the Kuiper Belt family proposed here. Within the former Koronis sub-family, the spread in free eccentricities is still larger than the spread in free inclinations, and appears to demand an anisotropic distribution of ejecta velocities. The corresponding widths for our candidate Kuiper Belt family are 0.15 and 0.03, respectively. That these widths are each $``$6 times greater than their Koronis sub-family counterparts probably reflects the greater ease with which orbital trajectories are altered in the more distant Kuiper Belt; orbital velocities are $``$4 times lower in the CKB than in the asteroid belt. The larger variation in $`e_{\mathrm{𝑓𝑟𝑒𝑒}}`$ compared to that in $`i_{\mathrm{𝑓𝑟𝑒𝑒}}`$ exhibited by both families presumably reflects a collision in which fragment ejecta velocities were greater in the (invariable) plane of the solar system than out of the plane. This, in turn, probably implicates a collision between two ancient bodies having greater relative velocity in the plane than out of the plane.

The mass of the projectile, $`M_{\mathrm{𝑝𝑟𝑜𝑗}}`$, that shattered the parent body could not have been much smaller than that of the parent, $`M_{\mathrm{𝑝𝑎𝑟}}`$. For the kinetic energy of the collision (measured in the frame moving at the center-of-mass velocity) to exceed the gravitational binding energy of the parent body, $`\mathrm{min}(M_{\mathrm{𝑝𝑟𝑜𝑗}}/M_{\mathrm{𝑝𝑎𝑟}})(v_{\mathrm{𝑒𝑠𝑐},\mathrm{𝑝𝑎𝑟}}/v_{\mathrm{𝑟𝑒𝑙}})^2`$, where $`v_{\mathrm{𝑒𝑠𝑐},\mathrm{𝑝𝑎𝑟}}`$ is the escape velocity from the parent body and $`v_{\mathrm{𝑟𝑒𝑙}}`$ is the relative velocity between projectile and target. Using the fragment sizes in Table 1 that are derived from an assumed albedo $`A=0.04`$, we estimate the minimum diameter of the putative parent body to be $``$$`800\sqrt{0.04/A}\mathrm{km}`$, for which $`\mathrm{min}(v_{\mathrm{𝑒𝑠𝑐},\mathrm{𝑝𝑎𝑟}})0.4\sqrt{0.04/A}\mathrm{km}\mathrm{s}^1`$. The relative velocity must be a fraction, $`f`$, of the local orbital velocity: $`v_{\mathrm{𝑟𝑒𝑙}}1(f/0.2)\mathrm{km}\mathrm{s}^1`$. Thus, $`\mathrm{min}(M_{\mathrm{𝑝𝑟𝑜𝑗}}/M_{\mathrm{𝑝𝑎𝑟}})0.2(0.2/f)^2(0.04/A)`$.

Surface reflectance spectra of members of a dynamical family can lend supporting evidence for a common physical ancestry. The Koronis sub-family exhibits relatively homogeneous S-type optical colors that set it apart from the more variegated field population (e.g., Gradie, Chapman, & Williams 1979). The parent body for the Koronis family is thought to be $``$$`120\mathrm{km}`$ in diameter (Bottke et al. 2001), too small to have undergone significant internal differentiation. More massive parent bodies may have differentiated; collisional fragments from such bodies may display greater variation in their spectral properties. Thus, finding a large spread in optical colors among members of a given candidate family does not necessarily eliminate the candidate from the running; widely different colors and/or spectra may instead offer a glimpse into the shattered, chemically zoned interior of a minor planet.

Unfortunately, color photometry exists for only 3 of our 32 candidate family members (1999 CM119, 1999 HU11, and 1999 CO153) in Table 1. Their B-V and B-R colors are plotted in Figure 2, along with the available colors of all other CKBOs \[we employ data from Trujillo & Brown (2002), selecting only those objects having $`43<a(\mathrm{AU})<47`$\]. The tendency for the 3 candidate family members to be redder than those of the general Classical population (Figure 2a) may imply homogeneity of the interior of the putative parent body. Alternatively, it may reflect whatever unknown process is responsible for generating the color-inclination correlation reported by Trujillo & Brown (2002). In an attempt to separate out this effect, we plot only data for those CKBOs having $`i_{\mathrm{𝑜𝑏𝑠}}<7^{}`$ in Figure 2b. There remains a tendency for our KBO family candidates to skirt only the lower envelope of the space spanned by all low-$`i`$ CKBOs, though clearly there are too few points and the uncertainties in individual points are too large to claim significant segregation. We urge physical studies of the candidate family members listed in Table 1 to further explore the possibility that we are indeed viewing the heretofore unseen interior of a massive, primitive planetesimal.

I am grateful to Marc W. Buie and Amy B. Jordan for inspiring this calculation during an observing run plagued by fog and hail. I thank also an anonymous referee for a thoughtful and thought-provoking report that significantly improved this paper. This work was supported in part by a Junior Faculty Research Grant awarded by the University of California at Berkeley.


text:
# Quantum Conductance in Semimetallic Bismuth Nanocontacts

## Abstract

Electronic transport properties of bismuth nanocontacts are analyzed by means of a low temperature scanning tunneling microscope. The subquantum steps observed in the conductance versus elongation curves give evidence of atomic rearrangements in the contact. The underlying quantum nature of the conductance reveals itself through peaks in the conductance histograms. The shape of the conductance curves at 77 K is well described by a simple gliding mechanism for the contact evolution during elongation. The strikingly different behaviour at 4 K suggests a charge carrier transition from light to heavy ones as the contact cross section becomes sufficiently small.

In the last few years, several techniques have been extensively used to study the conductance of nanometer-scale contacts with variable cross sections . Most of the previous works have been focused on metallic contacts, for which the Fermi wavelength $`\lambda _F`$ is very close to the dimensions of the atom. The elongation of the contact produces discrete changes in the cross section (i.e. in the number of atoms) resulting in abrupt jumps in the experimental $`GZ`$ curves (conductance vs. elongation of the contact). These jumps are of the order of the conductance quantum $`G_0=2e^2/h`$ since the addition of one atom to the contact is accompanied by the opening of one or few conduction channels depending on the metal. This interplay between electronic and mechanical properties in metallic contacts  makes it difficult to observe “pure” quantum effects on the conductance.

On the other hand, semimetallic nanocontacts could provide a unique laboratory to investigate quantum conductance phenomena in three-dimensional wires. The high value of $`\lambda _F`$ in these systems requires large contact areas in order to have total transmissions of the order of a conductance quantum. As a consequence, it will be possible to distinguish between atomic rearrangements and quantum effects. However, earlier experimental results did not clarify the situation. The first experiments  on semimetallic Sb contacts showed several sharp jumps in the $`GZ`$ curves for conductances much lower than the quantum unit, proving that, indeed, many atoms are necessary in order to obtain a value of the conductance equal to $`G_0`$. Nevertheless the conductance quantization features were absent probably due to a strong backscattering  or to the large angular opening  of the contact. More recently, conductance experiments on Bi nanowires at 4K presented a very different behaviour showing, in some situations, long conductance plateaus close to integer multiples of $`G_0`$. After statistical averaging of many consecutive curves, the histogram of conductance values showed the existence of two peaks close to $`G_0`$ and $`2G_0`$. However, the shape of the histogram exhibited strong deviations from those reported in metallic contacts.

In order to shine some light on the physics of semimetallic contacts, we have performed scanning tunneling microscopy (STM) experiments on Bi nanocontacts at different temperatures. First, we analyze in detail the opening of the first conduction channel, which is a particularly difficult task in metallic contacts, where the transition from the tunneling regime into contact is governed by the well known jump-to-contact phenomenon. The peculiar electronic structure of bismuth enables a selection of the carrier type, depending on the dimensions of the contact and temperature. This effect will be studied at 77K and 4.2K for contacts with typical conductances in the range $`0<G<2G_0`$. The coexistence of conductance quantization effects and atomic rearrangements will be studied through conductance histograms of the $`GZ`$ curves. Due to the fact that even for low conductance values many atoms will be forming the contact, the evolution of the $`GZ`$ curves will be analyzed in terms of a simple gliding mechanism that describes the variation of the contact shape during the elongation process.

The experiments have been performed using a high stability scanning tunneling microscope (STM) which can be operated over the range 2 K - 300 K inside a 4He cryostat. We have used high-purity (99.999%) polycrystalline bismuth, cut from a rod, for both STM tip and sample. Prior to the measurements the surfaces were cut and scratched. At a given location on the sample, the tip is crashed forcibly and repeatedly into the substrate, causing extended plastic deformations that lead to a connective neck between tip and sample. Different neck geometries can be obtained by varying the indentation depth . In Fig. 1 (Fig. 2 ) we show typical $`GZ`$ curves corresponding to a series of consecutive elongation cycles measured at 77K ( 4K ). Continuous cycles of elongation-contraction of a neck lead to situations in which the $`GZ`$ curve presents clear steps over the entire conductance range (typically from 0 to 10 G0). Although the curves are not reproducible they show a very similar shape. Typically, the length of the steps varies between 0.2 and 1 nm, depending on the initial conditions of the elongation-contraction cycle (i.e., the neck). These values are close to the lattice parameter of Bi (0.475nm). At 77 K the conductance steps have a small positive slope without specific structure. However, at 4K the first plateaus with $`GG_0`$ usually show smooth curved shapes which resemble those observed for metallic Al contacts . Moreover, the conductance discontinuities $`\mathrm{\Delta }G`$ are closer to $`G_0`$, suggesting that at 4K bismuth exhibits a metallic-like behaviour, with a smaller $`\lambda _F`$. Therefore, only a few atoms are involved in the contact. Distortions of the lattice parameter at the contact, together with the peculiar band structure of bismuth, are responsible of the curved plateaus obtained at 4 K.

As expected for a semimetallic contact, in the subquantum conductance regime the $`GZ`$ curves show clear steps at both temperatures (see Figs. 1 and 2). Except for the small jump discontinuities (caused by rearrangements of the atoms and thus reflecting mechanical contact), $`G`$ decreases exponentially as a function of $`Z`$ like in a vacuum tunneling process, but with an effective tunneling barrier that strongly depends on temperature (see Figs. 1b and 2b). While at 77K the effective barriers are typically of the order of 10meV, those obtained at 4K are close to 1eV. When the contact finally breaks, the conductance is usually so small that, due to limitations of our current experimental setup, we are not able to observe the transition from “contact” tunneling to actual vacuum tunneling. All these facts indicate that, even at these very small values of the conductance, tip and sample are still in mechanical contact. At 4K there are some situations (see Fig. 2(b), and the curve on the right hand side in Fig. 2(c)), in which the contact breaks before the subquantum conductance regime. In those cases, the conductance does not display subquantum steps and decreases exponentially as a function elongation, where the effective tunneling barrier $`\varphi 4`$eV is very close to the expected bismuth work function.

A detailed theoretical description of the electrical properties of the contact during elongation is a very difficult task. The effective mass of the charge carriers in Bi varies over a wide range from $`m^{}0.002m_0`$ for light electrons to $`m^{}m_0`$ for heavy electrons ($`m_0`$ is the electron mass) and depends on the crystal orientation. Moreover, the peculiar electronic properties of bismuth (semimetallic character, band filling, Fermi surface, …) may also vary with lattice deformations due to the extreme uniaxial strains occurring in small contacts during the elongation-compression processes. The apparent metallic behaviour at 4K can be qualitatively understood in terms of quantum finite size effects : as the contact radius $`R`$ decreases the quantized transversal momentum $`k_n1/R`$ increases. This leads to a shift in the subband energy position , $`E_n\mathrm{}^2/(2m^{})k_n^2`$ (larger for light carriers) and the gap between light and heavy electrons increases. At low temperatures, there is a transition from light carriers at large cross sections to predominantly heavy electron transport for small cross sections where the light electron states become empty. This qualitative discussion is fully consistent with the analysis of the tunneling regime in mechanical contact ($`G<G_0`$): The effective barriers in mechanical contact are related to the electron confinement in the constriction. In a free electron picture, the longitudinal momentum of each transversal mode or channel $`n`$ is given by $`k_{z_n}=\sqrt{K_F^2k_n^2}`$ (with $`K_F=2\pi /\lambda _F`$) such that to each channel $`n`$ corresponds a 1D effective potential barrier in the longitudinal direction . The modes with $`k_n>K_F=(2\pi /\lambda _F)`$ are evanescent or tunneling modes. When $`R\lambda _F`$ there is no room even for a single channel and the conductance is dominated by electron tunneling through the first channel with $`k_{z_1}ik_ni\pi /R`$. Assuming that the constriction length increases linearly with the elongation $`Z`$ we will roughly have $`G\mathrm{exp}(\frac{\pi }{R}Z)`$. Therefore, in the tunneling regime (before the contact breaks), the slope of $`\mathrm{ln}G`$ as a function of $`Z`$ gives the order of magnitude of the lateral size of the constriction. Our results then suggest, that, at high temperatures, subquantum steps are related to structural changes of constrictions involving hundreds of atoms ($`R`$ is a few nanometers). Since the Fermi wavelength is larger or of the order of $`R`$, it would also be of the order of a few nanometers. In contrast, at low temperatures subquantum conductance steps can be associated with changes in constrictions consisting of a few atoms ($`R`$ is now a few Å). This would imply $`\lambda _F`$ to be one order of magnitude smaller than before. As mentioned above, this qualitative discussion is fully consistent with the behaviour of the $`GZ`$ curves in the range of a few conductance quanta.

The observation of subquantum steps in the $`GZ`$ curves provides evidence of atomic rearrangements in the contact. Since these curves are not exactly reproducible from cycle to cycle, conductance quantization effects are difficult to observe from a single elongation curve. After the formation of a given neck, a histogram is obtained from 200 consecutive $`GZ`$ curves. In fig. 3(a) we present a set of histograms obtained from different initial necks. As it can be seen, the histograms present peaks close to integer multiples of $`G_0`$ which, in contrast with metallic contacts, can be unambiguously associated with conductance quantization effects. In general, we observe that the histogram shape (peak height, number of peaks, …) changes from neck to neck . This behaviour is clearly different than that observed in some metals (e.g. gold nanocontacts ). After averaging over 10 different necks, the resulting histogram (uppermost in Fig. 3a) still shows clear peaks at 1 and 2 conductance quanta which now resemble those of metallic contacts. However, in semimetallic contacts a variation of a few quanta in the conductance requires cross sections much larger than those in metallic ones. Therefore, the mechanical evolution in these systems is expected to be rather different than in metals. For cross sections larger than a few nanometers, yielding occurs as a consequence of collective shear events . A qualitative picture of the mechanical evolution can be obtained from a simple model system. We shall restrict ourselves to the tensile loading regime during the elongation of the neck. Suppose that the initial constriction is represented by a cylinder of radius $`R_0`$ and consider the deformation process sketched in the inset in Fig. 4(b). The contact is assumed to glide along a symmetry plane that forms an angle $`\theta `$ with respect to the cylinder axis. The maximum elongation $`Z_{max}`$ before the contact breaks down, is given by $`Z_{max}=2R_0/\mathrm{tan}(\theta )`$. For a finite elongation $`\stackrel{~}{Z}=Z/Z_{max}`$, the contact area $`A`$ shrinks according to (see inset Fig. 4(b)):

$$A(\stackrel{~}{Z})=2R_0^2\left(\text{acos}(\stackrel{~}{Z})\stackrel{~}{Z}\sqrt{1\stackrel{~}{Z}^2}\right).$$
(1)
Interestingly, the contact loses its axial symmetry and the narrowest section of the contact resembles an elliptical cross section (Fig. 4(b)) with an axis aspect ratio $`\eta (\stackrel{~}{Z})\pi R_0^2(1\stackrel{~}{Z})^2/A`$. An estimation of the transmission coefficients $`T_{nm}`$ as a function contact elongation can be obtained from a simple model based on the saddle-point-contact model . The geometry of our contact resembles the wide-narrow-wide geometry with ‘elliptical’ cross section discussed in Ref. . Therefore, $`T_{nm}`$ can be approximated by:

$$T_{nm}=\left[1+\mathrm{exp}\left(\frac{2\pi (\sqrt{\pi A}ϵ_{nm})}{\mathrm{tan}\alpha }\right)\right]^1$$
(2)
where $`ϵ_{nm}=(n+1/2)/\sqrt{\eta }+\sqrt{\eta }(m+1/2)`$ (all the lengths are in units of $`\lambda _F`$), and $`\mathrm{tan}\alpha =\sqrt{(A(0)A(\stackrel{~}{Z}))/\pi }/Z_{max}`$ defines the effective “opening angle”  of the constriction. The conductance curves are obtained from the Landauer formula, $`G=_{nm}T_{nm}`$.

Typical $`GZ`$ curves obtained by this simple model for $`\theta =50^{}`$ are shown in Fig. 4 together with some experimental results at 77K. In order to mimic the stepped structure, the cross section is assumed to be constant within a lattice parameter of Bi. Experimentally the conductance changes from $`10G_0`$ to $`G_0`$ after an elongation of the contact of the order of 15 nm. Within our model this fixes the Fermi wavelength to be in the nanometer scale, in agreement with the analysis of the subquantum regime. For small $`Z`$, i.e. for small opening angles, there is almost no scattering (the $`T_{nm}`$’s are 0 or 1). Therefore, the conductance is simply given by the number of propagating channels in the initial contact $`N_0`$, i.e., those $`nm`$ channels having $`ϵ_{nm}<\sqrt{\pi A}`$. As the elongation proceeds an interesting interplay between the narrowing of the constriction and the change of degeneracy of transversal modes can be observed. For $`\stackrel{~}{Z}=0`$ the contact has cylindrical symmetry and the transversal levels exhibit the well known degeneracy of a parabolic confining potential. Depending on the initial contact cross section, the conductance is then given by $`1,3,6,10,\mathrm{}G_0`$. As $`\stackrel{~}{Z}`$ increases, the cross section decreases and, at the same time, the degree of anisotropy $`\eta (\stackrel{~}{Z})`$ also varies, removing the mentioned degeneracy. The conductance behaviour along the elongation path depends critically on the initial conditions of the contact. The theoretical histograms (Fig. 3(b)) strongly depend on the initial cross section and exhibit an overall shape that is very similar to the experimental one.

In summary, the above experiments provide direct experimental evidence for the presence of conductance quantization in three-dimensional nanocontacts, as well as an experimental exploration of the conductance regime prior to the opening of the first quantum channel. Our results for moderate-sized contacts are well described by a simple model of contact evolution during an elongation process that is based on a gliding mechanism. In agreement with magnetotransport measurements, the differences observed between 77K and 4K curves for $`GG_0`$ are consistent with a transition from light to heavy electrons at the Fermi level when the contact cross section becomes sufficiently small. Finally, we consider that our work opens a new way to perform experiments on quantum coherence phenomena in atomic scale contacts that up to now were doable with semiconducting heterostructures only.

This research was supported by the DGICyT (Spain) through grants BFM2000-1470-C02-02, PB98-0464 and PB97-0068.


text:
# Eccentricity Excitation and Apsidal Resonance Capture in the Planetary System Upsilon Andromedae

## 1 INTRODUCTION

The outer two known planets orbiting Upsilon Andromedae ($`\upsilon `$ And) possess remarkably large orbital eccentricities (Butler et al. 1999). Planet C has a mass of $`m_1=(1.83/\mathrm{sin}i_1)M_J`$ (where $`i_1`$ is the angle between the orbit pole of planet C and our line of sight, and $`M_J`$ is the mass of Jupiter), an orbital semi-major axis of $`a_1=0.805`$ AU, and an orbital eccentricity of $`e_1=0.253`$. For planet D, the corresponding quantities are $`m_2=(3.79/\mathrm{sin}i_2)M_J`$, $`a_2=2.48\mathrm{AU}`$, and $`e_2=0.308`$. These orbital parameters were kindly supplied by D. Fischer (2002, personal communication), and represent more up-to-date values than those employed by previous works. Studies of the dynamical stability of the system (Rivera & Lissauer 2000; Stepinski, Malhotra, & Black 2000; Lissauer & Rivera 2001; Chiang, Tabachnik, & Tremaine 2001) strongly suggest that the two planets execute approximately co-planar trajectories that are observed nearly edge-on, so that $`\mathrm{sin}i_1\mathrm{sin}i_20.5`$.

That $`m_2`$ likely exceeds $`m_1`$ while $`e_2>e_1`$ stands at odds with the idea that gravitational interactions between planets D and C excited both $`e_2`$ and $`e_1`$ to their present-day values. Close encounters between planets (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Ford, Havlickova, & Rasio 2001) and mean-motion resonance crossings between planets on divergent orbits (Chiang, Fischer, & Thommes 2002) tend to impart the greater eccentricity to the less massive body. Thus, we are led to the conclusion that an external agent—another planet, a star, or the circumstellar disk from which the planets formed—must have played a direct role in exciting the eccentricity of the most massive planet, D. Of these three candidates, the latter is the most promising, as we explain below.

An additional clue to the origin of $`\upsilon `$ And lies in the apsidal resonance that obtains if $`i_{21}`$, the mutual inclination between the orbits of planets C and D, is less than $``$$`20\mathrm{°}`$.111The mutual inclination, $`i_{21}`$, does not necessarily equal $`i_2i_1`$; see, e.g., Stepinski, Malhotra, & Black (2000). The observation today that the arguments of periastron of the two orbits, $`\omega _1`$ and $`\omega _2`$, differ by only $``$$`12\mathrm{°}`$ has been used by Chiang, Tabachnik, & Tremaine (2001, hereafter CTT) to argue that the orbits are indeed nearly co-planar.222The argument of periastron, $`\omega `$, as fitted by Doppler velocity measurements is referred to the plane of the sky (see, e.g., CTT). The longitude of periastron, $`\stackrel{~}{\omega }`$, may be referred to any fixed line in inertial space for the purposes of this paper. When $`i_{21}=0\mathrm{°}`$, we can take $`\stackrel{~}{\omega }=\omega `$. If $`i_{21}20\mathrm{°}`$, then $`\mathrm{\Delta }\omega \omega _2\omega _1`$ librates about $`0\mathrm{°}`$ with a semi-amplitude of $``$$`38\mathrm{°}`$. This apsidal libration and the eccentricity evolution of both planets are displayed in Figure 1.

This paper lays the groundwork for understanding the origin of the large eccentricities and the apsidal alignment exhibited by the orbits of planets C and D in $`\upsilon `$ And. Our main result is that the eccentricity of planet C and the locking of orbital apsides are both consequences of the slow growth of the eccentricity of planet D. The latter eccentricity, in turn, was driven by an external agent—plausibly a primordial circumstellar disk lying exterior to the orbit of planet D—that acted over timescales exceeding $`10^4\mathrm{yr}`$. We play our scenario out and explain the mechanics of apsidal resonance capture in §2. We also discuss in that section why the recently detected companion star to $`\upsilon `$ And is not likely to be dynamically relevant. A summary is provided in §3.

## 2 Secular Growth of $`e_1`$ and Apsidal Resonance Capture

### 2.1 Basic Model

We take planets C and D to occupy co-planar, nearly circular orbits at time $`t=0`$. The dynamics exhibited by the system is qualitatively the same over the range $`0\mathrm{°}i_{21}30\mathrm{°}`$ (CTT). We adopt current observed values for the orbital semi-major axes, and minimum values for the planetary masses ($`\mathrm{sin}i_1=\mathrm{sin}i_2=1`$).

We introduce an external force only on planet D that smoothly amplifies that planet’s eccentricity over a finite time interval. The equations governing the eccentricities and apsidal longitudes ($`\stackrel{~}{\omega }`$) of the orbits of planets C and D read as follows (cf. Murray & Dermott 1999):

$`\dot{h_1}`$ $``$ $`{\displaystyle \frac{d}{dt}}(e_1\mathrm{sin}\stackrel{~}{\omega }_1)=+A_{11}k_1+A_{12}k_2`$ (1)
$`\dot{k_1}`$ $``$ $`{\displaystyle \frac{d}{dt}}(e_1\mathrm{cos}\stackrel{~}{\omega }_1)=A_{11}h_1A_{12}h_2`$ (2)
$`\dot{h_2}`$ $``$ $`{\displaystyle \frac{d}{dt}}(e_2\mathrm{sin}\stackrel{~}{\omega }_2)=+A_{21}k_1+A_{22}k_2+Eh_2`$ (3)
$`\dot{k_2}`$ $``$ $`{\displaystyle \frac{d}{dt}}(e_2\mathrm{cos}\stackrel{~}{\omega }_2)=A_{21}h_1A_{22}h_2+Ek_2`$ (4)
where

$`A_{11}`$ $`=`$ $`+{\displaystyle \frac{n_1}{4}}{\displaystyle \frac{m_2}{m_{}+m_1}}\left({\displaystyle \frac{a_1}{a_2}}\right)^2b_{3/2}^{(1)}\left({\displaystyle \frac{a_1}{a_2}}\right)`$ (5)
$`A_{12}`$ $`=`$ $`{\displaystyle \frac{n_1}{4}}{\displaystyle \frac{m_2}{m_{}+m_1}}\left({\displaystyle \frac{a_1}{a_2}}\right)^2b_{3/2}^{(2)}\left({\displaystyle \frac{a_1}{a_2}}\right)`$ (6)
$`A_{21}`$ $`=`$ $`{\displaystyle \frac{n_2}{4}}{\displaystyle \frac{m_1}{m_{}+m_2}}\left({\displaystyle \frac{a_1}{a_2}}\right)^1b_{3/2}^{(2)}\left({\displaystyle \frac{a_1}{a_2}}\right)`$ (7)
$`A_{22}`$ $`=`$ $`+{\displaystyle \frac{n_2}{4}}{\displaystyle \frac{m_1}{m_{}+m_2}}\left({\displaystyle \frac{a_1}{a_2}}\right)^1b_{3/2}^{(1)}\left({\displaystyle \frac{a_1}{a_2}}\right)`$ (8)
$`E`$ $`=`$ $`\mathrm{exp}(t/\tau _E)/\tau _e.`$ (9)
Here $`n_j`$ is the mean motion of planet $`j`$, $`m_{}=1.3M_{}`$ is the central stellar mass, $`b_{3/2}^{(1)}`$ and $`b_{3/2}^{(2)}`$ are the usual Laplace coefficients, and $`1/E`$ is the timescale over which exponential growth of $`e_2`$ occurs, where $`\tau _E`$ and $`\tau _e`$ are fixed time constants.

We emphasize that only the eccentricity of planet D is directly amplified by our external force. Our equations are intended to represent the following physical picture, staged just after the formation epoch of the planets. Interior to and between the orbits of planets C and D, circumstellar disk gas is absent. Only exterior to the orbit of planet D lies disk gas whose resonant interaction with that planet excites the planet’s eccentricity. Our scenario is consonant with most current thinking on planet-disk interactions. Gas between the two planets is plausibly driven out of this region by planet-induced torques, as numerical simulations by Kley (2000) and Bryden et al. (2000) find. Such gas either tunnels past both planets into the innermost and outermost disks, or is accreted by the planets (Bryden et al. 2000). Viscous accretion may drain away the disk that lies interior to both planets onto the central star, as has been invoked by Snellgrove, Papaloizou, & Nelson (2001). The disk outside $``$3 AU may be prevented from accreting inwards by shepherding torques exerted by the outermost planet. Goldreich & Sari (2002) have discovered that disk gas can excite a planet’s eccentricity; eccentricity-amplifying Lindblad resonances can defeat eccentricity-damping co-rotation resonances provided the planet’s initial eccentricity exceeds a critical threshold. While the nonlinear outcome of their finite amplitude instability remains to be worked out, it may be that a planet’s eccentricity ceases to grow when the planet’s orbit begins to overrun the gap edge, i.e., when $`e0.3`$. For a dimensionless disk viscosity of $`\alpha =10^4`$ and a disk mass of $`40M_J`$, we compute an eccentricity amplification timescale, $`e_2/\dot{e_2}`$, of $`7\times 10^4\mathrm{yr}`$, and a radial migration timescale, $`|a_2/\dot{a_2}|`$, of $`4\times 10^6\mathrm{yr}`$ \[Goldreich & Sari 2002; their equations (18) and (19)\]. Since the latter can exceed the former, we fix the semi-major axes of both planets in our analysis.333We have verified by direct numerical simulations that relaxing this assumption changes none of our conclusions.

We are aware that certain details underlying our scenario are not well understood. Gas between planets C and D might be shepherded rather than driven out (Bryden et al. 2000), in which case its effects on D and C would need to be accounted for. Numerical simulations by Papaloizou, Nelson, & Masset (2000) of interactions between Jupiter-mass planets and their parent disks result in eccentricity damping and have yet to be reconciled with the predictions of Goldreich & Sari (2002). Uncertainties in the plausibility of our scenario notwithstanding, the key features of equations (1)–(4) are that (a) the eccentricity evolution of planet D is directly driven by an external influence, and (b) the eccentricity and apsidal evolution of planet C are not, i.e., planet C feels only the gravitational potential of planet D. These minimal requirements are sufficient to understand the origin of the eccentricity of planet C and the origin of the curious apsidal alignment of the orbits of C and D, as we now demonstrate.

If $`E=0`$, equations (1)–(4) yield the classical Laplace-Lagrange (L-L) solution. One solution, fitted with initial conditions appropriate to observations of $`\upsilon `$ And today, is embodied in Figure 1. A variety of L-L solutions, differing only in initial eccentricities and apsidal longitudes, are mapped in the space of $`\mathrm{\Delta }\stackrel{~}{\omega }=\stackrel{~}{\omega }_2\stackrel{~}{\omega }_1`$ and $`e_2/e_1`$ in Figure 2a. The dashed separatrix divides circulating solutions from librating solutions. Librating solutions all have $`|\mathrm{\Delta }\stackrel{~}{\omega }|<90\mathrm{°}`$. Today planets C and D live on a contour of libration characterized by $`\mathrm{max}|\mathrm{\Delta }\stackrel{~}{\omega }|38\mathrm{°}`$.

If $`O(E)O(A_{jk})`$, or equivalently if $`1/E`$ is longer than apsidal precession timescales of order $`10^4\mathrm{yr}`$, then the system evolves adiabatically through a series of L-L solutions. Figure 2b follows the evolution of a system for which $`e_1(t=0)=e_2(t=0)=0.05`$, $`\mathrm{\Delta }\stackrel{~}{\omega }(t=0)=180\mathrm{°}`$, $`\tau _e=7.0\times 10^4\mathrm{yr}`$, and $`\tau _E=1.5\times 10^5\mathrm{yr}`$. The final time of the calculation is $`t_f=1\times 10^6\mathrm{yr}`$, at which time the eccentricity driving term, $`E`$, is effectively zero. Each point represents the position of the system in $`\mathrm{\Delta }\stackrel{~}{\omega }`$$`e_2/e_1`$ space at a particular time; the solution was obtained by numerically integrating equations (1)–(4) using a standard fourth-order Runge-Kutta routine (Press et al. 1992). A set of points highlighted by a symbol of a given type indicates the trajectory over a time interval of $`8400\mathrm{yr}`$. Comparison of Figure 2a with 2b indicates that the system morphs smoothly through a series of L-L contours. What is noteworthy is that the system is captured into libration and thereafter evolves towards smaller libration amplitude.

The mechanics underlying apsidal resonance capture and the damping of apsidal libration is as follows. The introduction of our eccentricity driving term amplifies $`e_2`$, i.e., it pushes the system to the right in Figure 2. Thus, a system initially occupying a contour of circulation is eventually pushed past the separatrix onto a contour of libration. It remains to explain the evolution towards smaller libration amplitude. A librating system executes an elliptical trajectory in $`\mathrm{\Delta }\stackrel{~}{\omega }`$$`e_2/e_1`$ space. While the system traces out the left-hand side of an elliptical contour, external amplification of $`e_2`$ pushes the system to the right, i.e., towards contours of smaller libration amplitude. While the system traces out the right-hand side of the ellipse, external amplification of $`e_2`$ also pushes the system to the right, but now towards contours of larger libration amplitude. Which half of the interaction dominates depends on how much time the system spends on each half of the ellipse. It is clear from Figure 1 that the system spends more time near the minimum value of $`e_2/e_1`$ than near the maximum value of $`e_2/e_1`$. The same behavior can be seen from Figure 2b; the points, plotted at equal time intervals, cluster more strongly near small $`e_2/e_1`$ than near large $`e_2/e_1`$. The reason for this is that the apsides of low eccentricity orbits tend to precess more quickly than those of high eccentricity orbits; for the same perturbative acceleration, $`\dot{\stackrel{~}{\omega }}\sqrt{1e^2}/e`$ by Gauss’s equation; and thus when $`e_1`$ is minimal ($`e_2/e_1`$ maximal), the orbital apsides of planet C precess more rapidly past the apsides of planet D than when $`e_1`$ is maximal ($`e_2/e_1`$ minimal). Variations in $`e_2`$ can be ignored in this analysis; they are characteristically smaller than variations in $`e_1`$ because planet D carries the lion’s share of the orbital angular momentum of the system (see, e.g., CTT and Figure 1). Thus, external driving of $`e_2`$ on the leftmost half of an elliptical contour dominates the driving of $`e_2`$ on the rightmost half, and the system spirals deeper into apsidal lock—all without the need for any explicit dissipation of energy!

The celestial mechanical considerations described above are general and we have verified by numerous integrations of equations (1)–(4) that apsidal resonance capture and the subsequent damping of apsidal libration occur under a wide variety of initial conditions. The input parameters employed for the trajectory shown in Figure 2b were chosen so that the end state of our calculation yielded system parameters that resemble those of $`\upsilon `$ And today. Figure 3 portrays the same model, displaying the eccentricity and apsidal evolution as explicit functions of time, and should be compared with Figure 1.

### 2.2 Refinements and Extensions

#### 2.2.1 Disk-induced Precession

In addition to exciting the eccentricity of planet D, an external disk would also cause the orbits of both planets to precess at rates different from those due solely to the planets’ mutual gravity. We may model disk-induced precession by adding extra terms to $`A_{11}`$ and $`A_{22}`$ that decay on timescales characterizing the dissipation of the disk (by, e.g., photoevaporation). This procedure can account for the disk’s temporary resonant and secular contributions to planetary precession rates. We have experimented by adding various transient terms to $`A_{11}`$ and $`A_{22}`$ and find that none of our conclusions changes qualitatively. Faced with extra, disk-driven differential precession, we may still reproduce the observed parameters of $`\upsilon `$ And by appropriately re-adjusting our input parameters (e.g., initial eccentricities and $`\tau _e`$).

#### 2.2.2 Stellar Companion

A stellar companion to $`\upsilon `$ And has recently been reported by 2MASS (Two Micron All-Sky Survey; Lowrance, Kirkpatrick, & Beichman 2002). The mass of the companion is $`m_E0.2M_{}`$ and the projected separation between it and the primary is 750 AU. Could this companion star be responsible for exciting the observed eccentricities of planets C and D?

The answer is almost certainly no. The usual Kozai mechanism for pumping planetary eccentricities (see, e.g., Holman, Touma, & Tremaine 1997) can only operate if (1) the plane of the planetary system is sufficiently inclined with respect to the plane of the stellar binary, and (2) the apsidal precession rate of the planet due to the stellar companion dominates other contributions to that precession rate. The latter condition fails to be satisfied, by a wide margin. We define $`e_E`$ and $`a_E`$ to be the (unknown) orbital eccentricity and semi-major axis, respectively, of the companion star. Then the precession period of planet D, as induced solely by the stellar companion, is of order $`2\pi n_2^1(m_{}/m_E)(a_E/a_2)^3(1e_E)^{3/2}2.5\times 10^8(0.2M_{}/m_E)(a_E/750\mathrm{AU})^3[(1e_E)/0.5]^{3/2}\mathrm{yr}`$ (Holman, Touma, & Tremaine 1997). This greatly exceeds the precession periods set by the mutual interaction between planets C and D, which are of order $`10^4\mathrm{yr}`$.

If the plane of the stellar binary has a small inclination with respect to the plane of the planetary system, we may treat the stellar companion as a third, massive body in our linearized equations of motion. It is not difficult to show that the companion star would then induce an eccentricity of order $`(a_2/a_E)e_E`$ in the orbit of planet D. For the same parameters employed above, this secularly forced eccentricity is of order $`10^3`$.

Sigurdsson (1992) has shown that the mean induced eccentricity and inclination of a planet orbiting a $`1.4M_{}`$ star exceed $`10^3`$ and $`10^2`$, respectively, if a $`0.7M_{}`$ companion approaches the primary to within a distance of $``$5 times the semi-major axis of the planet. These are mean induced orbital parameters, averaged over 25,000 different possible encounter geometries. For the stellar companion having $`a_E750\mathrm{AU}`$ to come within $`5a_213\mathrm{AU}`$ of the primary, $`e_E0.98`$. This is an underestimate of the eccentricity required for the companion star to be dynamically relevant because $`m_E0.2M_{}<0.7M_{}`$. In any case, we consider such large eccentricities to be improbable. Moreover, if such large eccentricities were to obtain, the consequent violent forcing of the planetary system over the orbital period of the star would threaten the stability of the system.

#### 2.2.3 Other Planetary Systems

We are aware of two other extrasolar two-planet systems whose present-day dynamics may be primarily secular: HD168443 (Marcy et al. 2001) and 47 UMa (Fischer et al. 2002). In the first system, the less massive, inner planet occupies the more eccentric orbit. This is consistent with the idea that both planets excited each other’s eccentricities by, e.g., divergent resonance crossings (Chiang et al. 2002). Alternatively, each planet’s eccentricity may have been excited by local disk material (Goldreich & Sari 2002), independently of the presence of the other planet. In HD168443, the separation in semi-major axis between the two planets is greater than that between planets C and D in $`\upsilon `$ And, increasing the likelihood that disk gas was confined between the two planets in the former system. In the absence of an innermost disk, ring confinement is necessary for the resonant excitation of the inner planet’s eccentricity, either by divergent resonance crossings or local disk-planet interaction. Both scenarios are consistent with the fact that HD168443 does not exhibit an apsidal lock; there is no reason to expect the ratio $`e_{outer}/e_{inner}`$ to be driven to large values.

In the case of 47 UMa, the eccentricity of the less massive, outer planet is constrained by observation to be less than 0.2 (Fischer et al. 2002). If said eccentricity is less than $`0.1`$, we find that the apsides of the two planets are locked about $`0\mathrm{°}`$, consistent with the findings of Fischer et al. (2002). Can we apply our theory for $`\upsilon `$ And to 47 UMa? The eccentricities in 47 UMa are so low that it is not clear how the system was driven by an external influence, if at all. Moreover, the phase space of possible trajectories in $`\mathrm{\Delta }\stackrel{~}{\omega }`$$`e_{outer}/e_{inner}`$ space is qualitatively different from the case of $`\upsilon `$ And, because most of the angular momentum of the system is carried by the inner planet. This causes an island of libration about $`\mathrm{\Delta }\stackrel{~}{\omega }=180\mathrm{°}`$ to appear at $`e_{outer}/e_{inner}1`$. An important feature of 47 UMa is the proximity of this system to a 5:2 mean-motion resonance; if future Doppler velocity data place the system in this resonance, then the apsides circulate (Fischer et al. 2002). Whether or not the system inhabits this mean-motion resonance, it seems that our theory for $`\upsilon `$ And cannot be applied to 47 UMa in a straightforward manner.

## 3 Summary

We have sketched an evolutionary scenario that explains the origin of the eccentricities and the apsidal alignment exhibited by the orbits of planets C and D in the extrasolar planetary system, Upsilon Andromedae. In our picture, the eccentricity of the outermost planet, D, grows smoothly via resonant interaction with a circumstellar disk of gas remaining from the formation of planets C and D. Provided the external, disk-driven amplification of the eccentricity of planet D occurs over a timescale that is long compared to apsidal precession timescales set by the planets’ mutual gravity, the orbits of planets C and D evolve adiabatically through a series of classical Laplace-Lagrange solutions. Increasing the ratio of the eccentricity of planet D to that of planet C eventually forces apsidal circulation to give way to apsidal libration. Continued pumping of this eccentricity ratio during the phase of the libration when this ratio is smallest drives the libration amplitude to zero. The resonantly excited eccentricity of planet D is secularly shared with that of planet C. In short, the substantial eccentricity of planet C and the apsidal resonance are both consequences of the growth of the eccentricity of planet D which, in turn, was driven by a third agent, most likely a remnant circumstellar disk lying exterior to the orbits of both planets. The companion M4.5V star to $`\upsilon `$ And recently detected by 2MASS is likely dynamically decoupled from the planetary system.

Our mechanism for capture into apsidal resonance and damping of apsidal libration is, to our knowledge, novel. Moreover, it enjoys complete independence from energy dissipation mechanisms that are normally required to damp libration amplitudes as in, e.g., the case of the Laplace resonance inhabited by the Galilean satellites (see, e.g., Murray & Dermott 1999). The gravitational interactions between planets C and D are, to excellent approximation, secular and therefore energy-conserving.

It is a pleasure to acknowledge Geoff Marcy and Debra Fischer for helpful and friendly exchanges and for providing orbital parameters, and Patrick Lowrance and Chas Beichman for generously sharing their most recent discovery.