text:
# Josephson vortex in a ratchet potential: Theory

## I Introduction

To extract useful work from random motion was a dream of mankind since the days when the Brownian motion was recognized. Unfortunately, the second law of thermodynamics forbids to extract energy “for free” from equilibrium thermal fluctuations (white noise), which was didactically demonstrated by R. Feynman in his Lectures. Nevertheless, one can extract useful work from non-equilibrium or time-correlated (colored) noise “not paying” for it, using so-called ratchets, i.e., systems with an asymmetric periodic potential. Recently there was a boost of activity in this field related to the experimental investigation of directed motion in biological systems, so-called Brownian motors which, e.g., move muscles or transport vesicles in a cell. In the latter case the probable mechanism of operation is the motion of kinesin molecules along the surface of microtubules, which can be mapped to the motion of Brownian particles along a one-dimensional ratchet potential with the period $`8.2\mathrm{nm}`$. The non-equilibrium energy is supplied by chemical reaction of splitting of adenosine triphosphate (ATP) which takes place close to the kinesin molecule.

In addition to the application of ratchets as noise rectifiers, it was suggested to use them for very efficient separation of small objects with different mobility, e.g., DNA molecules, viruses, etc.. The particle separation is based on so-called deterministic ratchets, where the particles move in a certain direction under the action of a deterministic force with zero time average. Moreover, changing the force profile one can reverse the direction of the particle motion. The classification and discussion of different types of ratchet systems can be found in Ref. .

In this paper we focus on Josephson ratchets which are of particular interest because (a) the directed motion results in a dc voltage according to the Josephson relation and (b) these systems can operate at very high frequencies up to about $`100\mathrm{GHz}`$. As a first example we mention the asymmetric dc SQUID where the equation of motion for the Josephson phase (difference of quantum mechanical phases) corresponds to the motion of an imaginary particle in a 2D ratchet potential. Such SQUID ratchets have been proposed  and studied experimentally. Another type of Josephson ratchet investigated recently is a 1D array of Josephson junctions with spatially modulated properties. A Josephson kink (vortex), which can move along the array, can be considered as quasi-particle in a 1D ratchet potential.

Here we propose a new class of Josephson ratchets which further develops the idea of a kink in a 1D array. The proposed system consists of a 1D long Josephson junction (LJJ) which may be bent in the $`ab`$-plane \[see Fig. 1(a)\] or have variable width $`w(x)`$. Here and below $`x`$ is a curvelinear coordinate along the junction. The fluxon (Josephson vortex) moving along the junction, from a mathematical point of view, is a topological soliton. It has its own mass, velocity, and other particle-like properties. We study the motion of a fluxon in LJJs in a ratchet potential which can be formed either by applying an external magnetic field and bending the junction properly or by modulating its width $`w(x)`$. To provide the required periodicity of the potential, the junction is topologically closed in a loop. Such a geometry is similar to the well known annular Josephson junction in which a fluxon moves in a sinusoidal potential created by a magnetic field. Using a more elaborated shape one can form an asymmetric potential with the possibility to control its amplitude by changing the amplitude of the external magnetic field. An alternative idea of magnetic field modulation using a specially shaped control line is proposed in Ref. . The directional motion of a fluxon can be detected by measuring the dc voltage across the junction which is, due to the Josephson relation, proportional to the average velocity of a fluxon.

Before discussing fluxon dynamics in a ratchet potential we would like to stress the difference between conventional ratchets and Josephson vortex ratchets. First, the fluxon, although it has a topological charge and in general behaves like a particle, is a nonlinear wave, i.e., it can change its shape rather strongly as well as emit electromagnetic waves. Second, the fluxon dynamics is usually studied in the underdamped limit which is opposite to the overdamped case which was considered for the majority of work on Brownian particles in a ratchet potential. Small damping may result in chaotic dynamics even in the deterministic case and, in the case of the fluxon, even without any potential. Thus, to have a well defined behavior of a fluxon ratchet, one has to work in the overdamped or in the weakly underdamped limit. If one uses conventional Nb-AlOx-Nb technology to fabricate the LJJ, this requirement means that the working temperature should be very close to $`T_c`$. As an alternative one can use junctions with intrinsically high damping such as SINIS LJJs or high-$`T_c`$ LJJs technology which allows to fabricate LJJs of required topology. Third, if we consider multi-particle dynamics, the strong repelling interaction between fluxons plays an important role and must be taken into account.

This paper is organized as follows. In section II we derive the equations for the dynamics of the Josephson phase in a bent LJJ of variable width $`w(x)`$ in the external magnetic field. We also discuss different kinds of fluxon ratchets, their advantages and drawbacks. The numerical simulation results are presented in III. Section IV concludes this work.

## II The Model

Here we derive the generalized perturbed sine-Gordon equation which takes into account the curvature of the LJJ in the $`ab`$-plane \[see Fig. 1(a)\], the uniform magnetic field $`\stackrel{}{𝐇}`$ applied in the plane of the junction (in $`b`$ direction) as well as the modulation of the LJJ width $`w(x)`$ along its length $`x`$.

We start from a discrete representation of the LJJ shown schematically in Fig. 1(b). The Kirchhoff equations for the Josephson phases in the cell and for the currents in one of the nodes are:

$`\varphi (x+dx)\varphi (x)={\displaystyle \frac{2\pi }{\mathrm{\Phi }_0}}\mathrm{\Phi }(x)={\displaystyle \frac{2\pi }{\mathrm{\Phi }_0}}\left[\mathrm{\Phi }_e(x)L(x)I_L(x)\right];`$ (1)
$`I_L(x)I_L(xdx)=I_e(x)I(x),`$ (2)
where $`\varphi (x)`$ is the Josephson phase at point $`x`$ of the junction, $`\mathrm{\Phi }(x)`$ and $`\mathrm{\Phi }_e(x)`$ are the total and the external magnetic flux applied to the cell, respectively, $`L(x)`$ is the inductance of the piece of the junction electrodes between $`x`$ and $`x+dx`$, $`I_L(x)`$ is the current in the electrodes, i.e. through the inductance $`L(x)`$, $`I_e(x)`$ is the externally applied bias current, and $`I(x,t)`$ is the current through the Josephson junction. The particular expression for $`I(x,t)`$ depends on the adopted JJ model and is introduced later.

Assuming that the interval $`dx`$ is infinitesimal, we can rewrite Eqs. (1) and (2) in a differential form using the following expressions:

$`I(x)`$ $`=`$ $`j(x)w(x)dx;`$ (3)
$`I_e(x)`$ $`=`$ $`j_e(x)w(x)dx;`$ (4)
$`L(x)`$ $`=`$ $`{\displaystyle \frac{\mu _0d^{}}{w(x)}}dx;`$ (5)
$`\mathrm{\Phi }_e(x)`$ $`=`$ $`\mu _0(\stackrel{}{𝐇}\stackrel{}{𝐧})\mathrm{\Lambda }dx=\mu _0H(x)\mathrm{\Lambda }dx,`$ (6)
where $`\mu _0d^{}`$ is the inductance of one square of the superconducting electrodes, $`d^{}2\lambda _L`$ is the effective magnetic thickness of the junction, $`\stackrel{}{𝐧}`$ is the unit vector normal to the plane of the junction cell as shown in Fig. 1(b), $`\mathrm{\Lambda }2\lambda _L`$ is the effective penetration depth of the magnetic field into the junction, and $`\lambda _L`$ is the London penetration depth of the superconducting electrode. We assume that the films are spatially uniform so that $`d^{}`$ and $`\mathrm{\Lambda }`$ are independent on $`x`$.

Substituting Eqs. (3)–(6) into Eqs. (1) and (2) we can rewrite the latter in a differential form as

$`{\displaystyle \frac{\varphi }{x}}={\displaystyle \frac{2\pi }{\mathrm{\Phi }_0}}\left[\mu _0H(x)\mathrm{\Lambda }{\displaystyle \frac{\mu _0d^{}}{w(x)}}I_L(x)\right];`$ (7)
$`{\displaystyle \frac{I_L(x)}{x}}=w(x)\left[j_e(x)j(x)\right].`$ (8)
Excluding $`I_L(x)`$ from Eqs. (7) and (8), we get the equation which describes the dynamics of the Josephson phase in the system:

$$w(x)\left[j_e(x)j(x)\right]=\frac{1}{\mu _0d^{}}\frac{d}{dx}\left\{w(x)\left[\mu _0\mathrm{\Lambda }H(x)\frac{\mathrm{\Phi }_0}{2\pi }\varphi _x\right]\right\},$$
(9)
Here and below, the subscripts $`t`$ and $`x`$, if any, denote the derivatives with respect to time $`t`$ and coordinate $`x`$, respectively. Note, that we didn’t include any particular model of JJ into our equation up to now, which is a definite advantage of this derivation procedure. In the case of the simple RSJ model, one should substitute $`j(x)`$ which is the sum of the supercurrent, normal (quasiparticle) current and displacement current densities:

$$j(x)=j_c\mathrm{sin}(\varphi )+\frac{\mathrm{\Phi }_0}{2\pi R}\varphi _t+C\frac{\mathrm{\Phi }_0}{2\pi }\varphi _{tt}$$
(10)
into Eq. (9). Here $`j_c`$, $`R`$ and $`C`$ are the critical current density, specific resistance and specific capacitance of the junction, accordingly. In this case Eq. (9) can be rewritten in a form which resembles the usual sine-Gordon equation:

$$\lambda _J^2\varphi _{xx}\omega _p^2\varphi _{tt}\mathrm{sin}(\varphi )=\omega _c^1\varphi _t\gamma (x)+QH_x(x)+\frac{w_x(x)}{w(x)}\left[QH(x)\lambda _J^2\varphi _x\right],$$
(11)
where $`\lambda _J=\sqrt{\mathrm{\Phi }_0/(2\pi \mu _0j_cd^{})}`$ is the Josephson penetration depth, $`\omega _p=\sqrt{2\pi j_c/(\mathrm{\Phi }_0C)}`$ is the Josephson plasma frequency, $`\omega _c=2\pi j_cR/\mathrm{\Phi }_0`$ is the so-called critical frequency, $`\gamma (x)=j_e(x)/j_c`$ is a normalized bias current density, and $`Q=2\pi \mu _0\mathrm{\Lambda }\lambda _J^2/\mathrm{\Phi }_0`$.

For theoretical investigation of the system we introduce standard normalized units, i.e., we normalize the coordinate to the Josephson penetration depth $`\lambda _J`$, and the time to the inverse plasma frequency $`\omega _p^1`$. After such simplifications, Eq. (11) can be rewritten as:

$$\varphi _{xx}\varphi _{tt}\mathrm{sin}(\varphi )=\alpha \varphi _t\gamma (x)+h_x(x)+\frac{w_x(x)}{w(x)}\left[h(x)\varphi _x\right],$$
(12)
with the damping coefficient $`\alpha =\omega _p/\omega _c1/\sqrt{\beta _c}`$, and the field $`h`$ normalized in a usual way as

$$h(x)=\frac{2H(x)}{H_{c1}}.$$
(13)
$`H_{c1}=\mathrm{\Phi }_0/(\pi \mu _0\mathrm{\Lambda }\lambda _J)`$ is the first critical field (penetration field) for a LJJ which is, in fact, equal to the field in the center of the fluxon. The normalized velocity is given in natural units of $`\overline{c}_0=\lambda _J\omega _p`$, where $`\overline{c}_0`$ is the so-called Swihart velocity. The normalized voltage $`V=\varphi _t`$ is given in units of $`\mathrm{\Phi }_0\omega _p/(2\pi )`$. From now on all quantities are given in normalized units.

In comparison with the usual perturbed sine-Gordon equation, Eq. (12) contains 3 additional terms. The term $`h_x(x)`$ describes the effect of the applied magnetic field when the junction is bent in the $`ab`$-plane. The second term $`\left[w_x(x)/w(x)\right]\varphi _x`$ comes from the width modulation. The last term $`\left[w_x(x)/w(x)\right]h(x)`$ describes the mixture of both and does appear only when both, field modulation due to curvature and width modulation are present.

It can be checked by a direct substitution into the Euler-Lagrange equation

$$\frac{d}{dt}\frac{}{\varphi _t}+\frac{d}{dx}\frac{}{\varphi _x}\frac{}{\varphi }=0,$$
(14)
that the Lagrangian density

$$=w(x)\left\{\frac{\varphi _t^2}{2}\frac{\left[\varphi _xh(x)\right]^2}{2}\left(1\mathrm{cos}\varphi \right)\right\},$$
(15)
results in the equation of motion (12) without the $`\alpha \varphi _t`$ and $`\gamma `$ terms which describe dissipation and external force and therefore are not included in the Lagrangian density. Of course, the Lagrangian density (15) can be obtained directly from Fig. 1(b) and the RSJ model. From Eq. (15) one can see that instead of the usual potential energy term $`\varphi _x^2/2`$ we now have $`w(x)\left[\varphi _xh(x)\right]^2/2`$, i.e., actually 3 terms. The first, $`w(x)\varphi _x^2/2`$, is the obvious generalization of the usual potential energy term to the case of variable width $`w(x)`$. The second term $`w(x)h(x)^2/2`$ is a constant term due to the applied field and is not related to the fluxon motion or other Josephson phase activity in the junction. In fact, there are no traces of this term in Eq. (12). The third term $`h(x)\varphi _xw(x)`$ represents the part of the potential energy density which we are interested in and which we are going to exploit to build a system with a ratchet potential.

One of the solutions of the sine-Gordon Eq. (12) with zero r.h.s. is a soliton (fluxon)

$$\varphi (x,x_0)=4\mathrm{arctan}\mathrm{exp}(xx_0),$$
(16)
with the center situated at point $`x_0`$. We consider non-relativistic motion i.e., $`dx_0/dt1`$. Further assuming that the fluxon profile (16) does not change much due to the r.h.s. of Eq. (12) which acts as a perturbation, we can get the explicit expression for the potential energy $`U(x_0)`$ as a function of fluxon coordinate $`x_0`$. For this purpose, we note that in expression (15) for the Lagrangian density the second and the third terms correspond to the potential energy density $`𝒰(x,x_0)`$ (with opposite sign). The potential energy $`U(x_0)`$ is obtained by substituting $`\varphi (x,x_0)`$ from Eq. (16) into $`𝒰(x,x_0)`$ and integrating over $`x`$. Thus, we get

$$U(x_0)=_{\mathrm{}}^+\mathrm{}\frac{4w(x)}{\mathrm{cosh}^2(xx_0)}\frac{2w(x)h(x)}{\mathrm{cosh}(xx_0)}dx.$$
(17)
The first term corresponds to the potential energy due to width modulation, the second due to shape, field and width.

The first possibility to form a ratchet potential is to apply no magnetic field ($`h=0`$) and to vary the width $`w(x)`$ of the junction. In this case the potential energy will be given by the first term of Eq. (17). Moreover, when the junction width $`w(x)`$ does not change much over the distance comparable with the fluxon size, the potential

$$U(x_0)8w(x_0)$$
(18)
just repeats the $`w(x_0)`$ profile. From Eq. (18) it follows that if the width changes as a saw-tooth, so does the potential energy, except for the vicinity of the saw-tooth’s infinite slope where one has to calculate the convolution according to Eq. (17). Here and below when mentioning saw-tooth profile we mean a saw-tooth with finite positive slope and infinite negative slope. In the case of the saw-tooth where the width changes from $`w_0`$ to $`w_0+\mathrm{\Delta }w`$ at the point $`x_0=0`$ (infinite slope), the potential energy in the vicinity of the point $`x_0=0`$ is $`U(x_0)=8w_0+4\mathrm{\Delta }w\left[1+\mathrm{tanh}(x_0)\right]`$. The practical implementation of such a ratchet would look like an annular LJJ with the outer edge having the shape of a circle and inner edge having the shape of one turn of a spiral. This geometry has the advantage that the corresponding potential can be made ideally saw-tooth-like except for the smearing due to the convolution in Eq. (17), which is a common feature of all fluxon based systems. In some sense, this system is an analog of the Josephson ratchets based on 1D arrays. Unfortunately both types of ratchets do not allow to control the potential height during experiment which can be considered as a disadvantage.

The second possibility is to keep $`w`$ constant, to apply a magnetic field and to bend the junction in the $`ab`$-plane. In this case the first term of Eq. (17) gives a constant and we have to consider only the second term. Again, if $`h(x)`$ changes slowly in comparison with the fluxon size, Eq. (17) is simplified to

$$U(x_0)=2\pi wh(x_0).$$
(19)
In the well known case of the ring shaped junction, the field $`h(x)=h_0\mathrm{cos}(\theta )=h_0\mathrm{cos}(x/R)`$ \[$`\theta `$ is the angle between $`\stackrel{}{𝐇}`$ and $`\stackrel{}{𝐧}`$, as shown in Fig. 1\], and therefore we get a symmetric potential. If we deform the ring properly, the potential can be made asymmetric as desired. Possible experimental shapes are shown in Fig. 2. The advantage of this kind of ratchets is that one can control the amplitude of the potential during experiment by varying the amplitude $`h_0`$ of the magnetic field. The possibility to tune the potential height allows to implement so-called flashing ratchets by applying either an ac magnetic field using a coil or just placing an rf antenna close to the junction so that the $`\stackrel{}{𝐇}`$ of the emitted electromagnetic wave will have non-zero component in the $`ab`$-plane.

Note, that the terms $`\gamma (x)`$ and $`h_x(x)`$ in Eq. (12) from a mathematical point of view play the same role. Therefore in the experiment the field $`h(x)`$ can be substituted by a properly chosen additional bias current $`\gamma _p(x)=h_x(x)`$, which has zero average in space. The inverse mapping is also valid, but the bias current with non-zero average maps to a non-periodic field (potential) with linearly growing component, which does not belong to the class of ratchets. In section III.A we, in fact, use $`\gamma 0`$ (these brackets denote spacial averaging), but only to test the asymmetry of the potential. When we demonstrate the real operation of the ratchet in the section III.B, we have $`\gamma =0`$.

## III Simulation results

In this section we study the fluxon dynamics in a Josephson ratchet of the second type ($`h0`$, $`w=\mathrm{const}`$) for a saw-tooth field profile $`h(x)`$. This ideal asymmetric profile is not only of academic interest because it should show the best figures of merit for fluxon ratchets but also can be quite closely reproduced in a real system \[see shape in Fig. 2(b) and corresponding $`h(x)`$ in Fig. 3(a)\]. All simulations were performed using an explicit numerical scheme for Eq. (12) using a LJJ of the normalized length $`\mathrm{}=20`$, with damping coefficient $`\alpha =0.2`$ (weakly underdamped limit). The numerical technique and simulation software are discussed in detail elsewhere.

### A Probing the asymmetry of the potential

First, we probe the asymmetry of the potential by calculating critical current $`\gamma _c`$ vs. potential height $`h_0`$ for the case of one trapped fluxon \[$`\gamma _{c1}^\pm (h_0)`$\] and for the case of no trapped fluxons \[$`\gamma _{c0}^\pm (h_0)`$\]. The superscripts ‘$`+`$’ or ‘$``$’ correspond to opposite directions of applied bias current.

On the basis of the model derived in the previous section we can understand how these dependences should look like for arbitrary magnetic field profile $`h(x)`$. In general, due to the loop-like geometry, the left and right tails of a fluxon can interact. We assume that the junction is long enough and this interaction in negligible. This situation is equivalent to the long periodic system where the fluxons are separated by a large distance $`\mathrm{}\lambda _J`$. Furthermore, we represent the field $`h(x)`$ as $`h_0(x)`$ and map it to the equivalent additional bias current $`\gamma _p(x)=h_0_x(x)`$. An example of $`h(x)`$ and $`\gamma _p(x)`$ derived from the geometry in Fig. 2(b) is shown in Fig. 3(a).

If there is a fluxon in the junction, $`\gamma (x)`$ and $`\gamma _p(x)`$ translate into a driving force $`F_\gamma (x_0)`$ and a potential (pinning) force $`F_p(x_0)`$, respectively, acting on the fluxon

$$F_\gamma (x_0)=_{\mathrm{}}^+\mathrm{}\frac{2\gamma (x)w}{\mathrm{cosh}(xx_0)}𝑑x=2\pi \gamma w,$$
(20)
$$F_p(x_0)=_{\mathrm{}}^+\mathrm{}\frac{2\gamma _p(x)w}{\mathrm{cosh}(xx_0)}𝑑x=h_0f(x_0)w.$$
(21)
The corresponding potential $`U(x_0)`$ and force $`f(x_0)`$ are shown in Fig. 3(b). For the sake of simplicity we suppose that $`\gamma `$ does not depend on $`x`$, but our discussion can be easily generalized to the case when the r.h.s. of Eq. (20) is equal to $`2\pi \mathrm{\Gamma }(x_0)w`$. The pinning force (21) can be also obtained directly from expression (17) for the potential energy. When we increase $`\gamma `$ the fluxon is pinned while these two forces can compensate each other, i.e., $`F_\gamma (x_0)+F_p(x_0)=0`$. The depinning happens for

$$\gamma =\gamma _{c1}=h_0f(x_1)/2\pi ,$$
(22)
where $`x_1`$ is the coordinate at which $`f(x_1)`$ has a minimum \[see Fig. 3(b)\]. We assume that $`\gamma >0`$ and $`f(x_1)<0`$. Thus from Eq. (22) we see that $`\gamma _{c1}(h_0)`$ looks like a straight line which starts from the origin.

On the other hand, regardless of the presence of a fluxon, the Josephson phase evolves under the action of the net current $`\gamma +\gamma _p`$. From Fig. 3(a) it is clear that at point $`x_2`$, where $`\gamma _p(x)`$ is maximum, this sum can exceed 1 for some value of $`\gamma =\gamma _c`$ and the junction switches to the resistive state. Using the $`\gamma _p=h_0_x(x)`$ the condition $`\gamma _c+\gamma _p(x_2)=1`$ gives

$$\gamma _c(h_0)=1+h_0_x(x_2)=1h_0\left|_x(x_2)\right|,$$
(23)
i.e., also a straight line, but with the negative slope.

Summarizing, $`\gamma _{c1}(h_0)`$ has two branches (22) and (23), and follows the one with the lowest critical current for given $`h_0`$. For fields $`0<h_0<h_{}`$, $`\gamma _{c1}(h_0)`$ is given by Eq. (22), and for $`h_0>h_{}`$ by Eq. (23), where the field

$$h_{}=\frac{1}{f(x_1)/2\pi +_x(x_2)},$$
(24)
is the field where these two dependences intersect. The dependence $`\gamma _{c0}(h_0)`$ has only one branch (23) since the fluxon depinning mechanism is absent.

When one applies the current $`\gamma `$ in the opposite direction, the dependences qualitatively look the same but the particular values of slopes can be different in the case of an asymmetric potential, if $`\gamma _p(x_2)\gamma _p(x_4)`$ and $`f(x_1)f(x_3)`$ (see Fig. 3). Thus, the measurement of $`\gamma _{c0}^\pm (h_0)`$ and $`\gamma _{c1}^\pm (h_0)`$ gives direct information about the asymmetry of the field and of the potential, respectively. Note that an asymmetric potential may have $`\gamma _p(x_2)=\gamma _p(x_4)`$ and/or $`f(x_1)=f(x_3)`$, i.e., it does not reveal its asymmetry in the measurements of $`\gamma _{c0}^\pm (h_0)`$ and $`\gamma _{c1}^\pm (h_0)`$. Inversely, the asymmetry of $`\gamma _{c0}^\pm (h_0)`$ and $`\gamma _{c1}^\pm (h_0)`$ is a clear indication of asymmetry in the system.

We stress that the first branch of $`\gamma _{c0}^\pm (h_0)`$ and the second branch of $`\gamma _{c1}^\pm (h_0)`$ show the asymmetry of the field $`h(x)`$, while the first branch of $`\gamma _{c1}^\pm (h_0)`$ shows the asymmetry of potential $`U(x_0)`$, i.e., already after convolution (17) with $`\varphi _x`$. While $`h_x(x_0)`$ may have appreciable asymmetry, $`U(x_0)`$ may be almost symmetric.

We should note that the above analysis is valid for an arbitrary field profile $`h(x)`$ as long as $`h_x(x)`$ can be considered as a perturbation in the r.h.s. of Eq. (12). However, for the case of the saw-tooth potential this is not the case, since $`h_x(x)`$ has a $`\delta (x)`$-like behavior at the point of the junction where the saw-tooth has infinite slope. Therefore our analysis is valid only qualitatively. In particular the dependences $`\gamma _c^\pm (h_0)`$ still have linear slopes but the values of the slopes can not be calculated using simple integrals as shown above because the fluxon shape will differ from (16) quite considerably, as we saw in numerical simulations. Therefore we present the results of direct numerical simulation of Eq. (12).

The simulated $`\gamma _{c0}^\pm (h_0)`$ and $`\gamma _{c1}^\pm (h_0)`$ for an ideal saw-tooth field profile $`h(x)`$ are presented in Fig. 4. The magnetic field is given in terms of the force amplitude $`h_02\pi /\mathrm{}`$. If the system is closed in a loop with a normalized circumference $`\mathrm{}`$, the $`h(x)`$ can be written in the form $`h_0(2\pi x/\mathrm{})`$ where $``$ is geometry dependent $`2\pi `$ periodic function normalized to 1. Therefore the term $`h_x(x)`$ in Eq. (12) is

$$h_x(x)=h_0\frac{2\pi }{\mathrm{}}^{}\left(\frac{2\pi }{\mathrm{}}x\right),$$
(25)
where $`^{}(\xi )=d(\xi )/d\xi `$ is a periodic function. This means that the potential force scales inversely proportional to the length of the junction $`\mathrm{}`$. So, to get rid of the length dependence and make our results valid for any length $`\mathrm{}1`$, we present our results as a function of $`h_02\pi /\mathrm{}`$.

From Fig. 4 one can see that the depinning current $`\gamma _{c1}(0)=0`$ and then grows linearly with field, as expected from the theory, up to some field $`h_{}2\pi /\mathrm{}1.8`$. After that it decreases linearly with field, also according to our prediction, and, finally, exhibits multiple branches corresponding to the multi-fluxon states.

As it should be in the case of a saw-tooth potential, the slopes of $`\gamma _{c0}^\pm (h_0)`$ and $`\gamma _{c1}^\pm (h_0)`$ are not symmetric, and the ratio of the slopes is about $`4`$. Ideally, if the fluxon were a real particle we could expect an infinite force necessary to push the fluxon out of the well in the direction against the infinite slope and, therefore, an infinite ratio of slopes on $`\gamma _{c1}^\pm (h_0)`$ curves. In practice, due to the convolution (17), we have a finite force and a finite ratio of slopes on the $`\gamma _{c1}^\pm (h_0)`$ dependences. The finite ratio of slopes on $`\gamma _{c0}^\pm (h_0)`$ as well as of the second slopes $`\gamma _{c1}^\pm (h_0>h_{})`$ is observed because perturbation theory does not apply as mentioned above.

The next step is the understanding of the dynamics of the fluxon in an asymmetric potential. One of the simplest observations which can be made numerically as well as experimentally is the examination of the fluxon trapping current $`\gamma _{\mathrm{tr}}`$, i.e., the minimum current at which a fluxon still moves along the system, not being trapped by the potential. Obviously, in the underdamped system the trapping of the fluxon by the potential will not take place while the kinetic energy of the fluxon exceeds the height of the potential, and one should not see any difference in the fluxon trapping currents $`\gamma _{\mathrm{tr}}^\pm (h_0)`$ for opposite bias current directions. In the strongly overdamped case the fluxon dissipates energy so fast that to move further under the action of the driving force, the driving force should always overcome the maximum value of the potential (trapping) force. In an asymmetric potential the maximum force created by the potential is different for opposite directions of fluxon motion and we should expect a difference in the fluxon trapping current. The simulated $`\gamma _{\mathrm{tr}}^\pm (h_0)`$ curves for a saw-tooth potential are shown in Fig. 5. Even in the slightly underdamped case $`\alpha =0.2`$ the dependences for positive and negative direction of fluxon motion (bias current) differ quite considerably. We also note that $`\gamma _{\mathrm{tr}}^+(h_0)`$ almost coincides with $`\gamma _{c1}^+(h_0)`$ (see Fig. 4) for $`h_0<2`$ which means that for positive direction of bias there is no hysteresis on the $`I`$$`V`$ characteristic when the system switches from zero voltage state to the state with moving fluxon and back. Instead, for the negative bias the hysteresis is present as can be seen in Fig. 6.

### B Motion due to monochromatic force

In the following we investigate the motion of a fluxon in a system with the saw-tooth field profile $`h(x)`$ and $`w=\mathrm{const}`$ under the action of a monochromatic ac bias current $`\gamma =\gamma _{\mathrm{ac}}\mathrm{sin}(\omega t)`$ in Eq. (12). If as a result of an ac drive the fluxon starts to move around the junction, we can estimate the characteristic frequency of this process. It is clear that the maximum velocity of the fluxon is equal to the Swihart velocity (1 in normalized units) therefore the maximum revolution frequency is $`\omega _{\mathrm{max}}=2\pi /\mathrm{}0.314`$.

In the limit of low frequency (quasi-static) drive $`\omega \omega _{\mathrm{max}}`$ we can calculate the net average velocity of the fluxon, just integrating the current-velocity characteristics $`u(\gamma )`$ (see Fig. 6) at a given value of field $`h_0`$ as:

$$u(\gamma _{\mathrm{ac}})=\frac{1}{2\pi }_0^{2\pi }u(\gamma _{\mathrm{ac}}\mathrm{sin}(\tau ))𝑑\tau ,$$
(26)
where the integration path should follow the proper hysteretic branches of the $`u(\gamma )`$ dependence.

In the case when $`\omega `$ is comparable with the maximum frequency of the fluxon rotation around the junction $`\omega _{\mathrm{max}}`$, we perform direct numerical simulations. The dependences $`u(\gamma _{\mathrm{ac}})`$ calculated numerically for different values of $`\omega `$ are shown in Fig. 7(a)–(c) for $`h_02\pi /\mathrm{}=0.5`$, 1.0 and 2.0, respectively. The upper curves marked as “qs” (quasi-static) are calculated using Eq. (26). All other curves for $`\omega =0.01`$–0.1 were obtained by means of direct numerical simulation of Eq. (12). Such values of frequencies were chosen taking into account that they should be less but comparable with $`\omega _{\mathrm{max}}`$. As we saw in simulations, at $`\omega 0.2`$ \[in the case of Fig. 7(a) already for $`\omega =0.1`$\] the rectification is suppressed almost completely.

From Fig. 7 one can see that the dependence $`u(\gamma _{\mathrm{ac}})`$ for quasi-static case is piecewise with jumps at the values of bias current where switching between different hysteretic branches occurs. For example, let us focus on $`h_02\pi /\mathrm{}=1.0`$ \[see Fig. 7(b)\]. In this case the switch between different hysteretic branches occurs at $`\gamma `$ which is equal to 0.522, 0.696, 0.916 (c.f. Figs. 4, 6 and 7). Several characteristic regions can be distinguished. First, in the region of small values of $`\gamma _{\mathrm{ac}}`$, where $`u=0`$, the amplitude of the ac current is not big enough to exceed the pinning current and the fluxon is localized in the well. At higher ac bias up to $`\gamma _{\mathrm{ac}}=0.522`$ the positive depinning current is exceeded and the system switches to the fluxon step, i.e., the fluxon escapes from the well and starts to move around the junction. Note that during a period of ac drive the fluxon can only make integer number of turns and at the end of the period will be again localized in the well. This results in the quantization of the average velocity $`u`$, as can be seen as steps on all curves except for the quasi-static one. The step number corresponds to the number of turns done by the fluxon during one period of ac drive. The voltage rectifier based on this principle will give a quantized voltage $`V_n=n\mathrm{\Phi }_0\omega /2\pi `$ (in physical units) which is defined by fundamental constants and the applied frequency. The corresponding velocity quanta in Fig. 7 are given by $`u_n=n\omega \mathrm{}/(2\pi )`$. This gives $`u_10.0318`$ for $`\mathrm{}=20`$ and $`\omega =0.01`$. Note that the voltage quantum increases with $`\omega `$ and for $`\omega =0.01`$ five quanta (turns per ac cycle) give the same voltage as one quantum for $`\omega =0.05`$. For the quasi-static case the voltage quantum is infinitesimal so that the quasi-static curves look smooth.

A further increase of the ac bias amplitude $`0.522<\gamma _{\mathrm{ac}}<0.696`$ results in a decrease of the dc voltage because during the negative semi-period of the ac drive the fluxon starts moving in the opposite direction so that the efficiency of the fluxon ratchet drops.

The big peak at $`0.696<\gamma _{\mathrm{ac}}<0.916`$ in Fig. 7(b) is related to the switching of the junction to the resistive state during the positive semi-period while during the negative semi-period the junction stays in the zero-voltage state or on the fluxon step. When the amplitude of ac drive gets big enough to switch the system to the resistive state also during negative semi-periods, the dc voltage drops again and in the limit of very strong ac drive approaches zero. A large amplitude of ac bias $`\gamma _{\mathrm{ac}}1`$ implies that the system spends almost the whole period in the positive or negative resistive state and only a tiny fraction of the period on the asymmetric part of the current-velocity characteristic at low currents, so that the resulting average velocity is close to zero. In addition, such a strong ac drive results in chaotic dynamics, as can be seen in Fig. 7 for large values of $`\gamma _{\mathrm{ac}}`$.

Figure 7 also shows that the rectification effect decreases when the driving frequency increases and approaches the maximum fluxon rotation frequency $`\omega _{\mathrm{max}}`$.

We point out that the performance of the fluxon rectifier depends on the chosen potential amplitude $`h_0`$. For comparison in Fig. 7(a)–(c) we show the rectification characteristics for different values of magnetic field amplitude $`2\pi h_0/\mathrm{}=0.5`$, 1.0 and 2.0 (see Fig. 4). We see that the potential depth affects the amplitude of ac drive at which rectification appears, at which $`u(\gamma _{\mathrm{ac}})`$ has maximum and it affects the relative location of the two regions corresponding to the fluxon rectification regime and resistive rectification regime. We note, that the resistive rectification regime in principle gives a larger effect which is inversely proportional to the damping coefficient $`\alpha `$, i.e., it can be made rather big. However, this is a rather trivial effect which can be obtained using any non-symmetric $`I`$$`V`$ curve. The fluxon rectification regime, which is the main subject of this study, gives smaller but a quantized voltages. This can be an advantage for some applications.

## IV Conclusion

We proposed a new type of Josephson vortex ratchets, where the motion of a fluxon along a long Josephson junction closed in a loop can be considered as the motion of a quasi-particle in a ratchet potential. The derived model suggests several implementations of the fluxon ratchets, and, in particular, the one where the amplitude of the potential can be controlled during experiment. Since the fluxon is a soliton (non-linear wave) moving in the underdamped medium with asymmetric potential, we expect some non-trivial effects related to (a) wave properties of the fluxon, (b) underdamped dynamics and (c) interaction between fluxons. As a first step, we performed a numerical analysis of the effective potential seen by a fluxon, and checked the rectification of a monochromatic signal. We found that a fluxon rectifier produces a quantized voltage, with a quantum given (due to the Josephson relation) only by the fundamental constant $`\mathrm{\Phi }_0`$ and the driving frequency $`\omega `$. The experimental investigation of the proposed system based on (Nb-Al-AlOx)-Nb Josephson junction technology is in progress.

###### Acknowledgements.

We thank P. Caputo, A. Ustinov, M. Fistul and B. Malomed for fruitful discussions.


format: text

text:
# Orbit Fitting and Uncertainties for Kuiper Belt Objects

## 1 Introduction

The discovery of the Kuiper Belt (Jewitt and Luu, 1993) created a new vista upon the formation and early evolution of the Solar System. Study of these objects’ basic physical properties of course requires sufficiently well-determined orbital parameters to deduce a future position and the distance to the object. The orbital parameters are themselves of great interest, with the phase space of orbits threaded by an intricate web of resonant and long-lived chaotic zones (Malhotra, 1996). For both the practical goal of object retrieval and the higher goal of understanding the dynamical structure of the population, it is essential to not only have accurate orbits but also to quantify the uncertainty in the orbit.

There are centuries-old methods for determination of orbital parameters from observed positions. Present-day desktop computers execute these solutions trivially, and a linearized propagation of the positional uncertainties into the orbital elements is straightforward as well (Muinonen and Bowell, 1993). Workstations are even fast enough to bound the uncertainty region with a brute-force sampling of the 6-dimensional orbit space in many cases. Why then should we bother developing another orbit-fitting technique? KBOs pose a particular challenge because the objects are faint, so recovery observations are quite costly, requiring substantial investment of 2–4-meter telescope time for all but the brightest known objects. As a consequence, most known objects have been observed only a few times, leading to substantial degeneracies in the orbit fits. We need to know the error ellipse (in sky position and in orbital parameter space) even in the face of these degeneracies. In the best case a brute-force method will work on underdetermined orbits but leave us ignorant of the nature of the degeneracy. In the worst case the degeneracies will lead to long or effectively infinite computation times.

The Minor Planet Center (MPC) (http://cfa-www.harvard.edu/iau/mpc.html) provides orbital elements and predicted ephemerides for new objects. The approach of the MPC to short-arc objects is to select the simplest “sensible” orbit that fits the data—e.g. a circular orbit, if it fits the data and does not imply a close approach to Neptune; a 3:2 resonant orbit with Neptune if the circular orbit would imply a close approach; or, if the circular orbit is a poor fit, a “Väisälä” orbit, with the object near perihelion (Marsden, 1985, 1991), if this does not imply a close encounter. The MPC orbits in effect remove the degeneracies by assuming that the orbit is most likely to resemble those of known objects, and consequently predict positions that in most cases are quite close to the actual recovery positions. There are a few difficulties with this approach, however: first, focussing recovery efforts on the positions predicted by these favored orbits will inevitably bias the recovered population toward such orbits. This makes it difficult to quantify the portion of the population which may be in unusual orbits as they are preferentially lost. Second, uncertainties on positions and orbital elements are not currently available from the Minor Planet Center. Positional uncertainties are essential for proper planning and evaluation of recovery observations. Uncertainties on orbital elements are important for analyses of the dynamical characteristics of the population. For example, there was speculation that the 2:1 resonance with Neptune was empty or underpopulated, but 1997 SZ10 and 1996 TR66 were reclassified from 3:2 resonators to 2:1 resonators after further observations. Finally, it is sometimes necessary to calculate orbits “on the fly” during an observing run, and the MPC should not be expected to provide instantaneous analyses.

## 2 Methods

### 2.1 Motivation

Our approach to orbit-fitting for KBOs is motivated and guided by the following differences from the historically more common application to the asteroid population:

* Recovery observations are very expensive, as noted. This leads us to investigate, in §4, the minimal number of observations required to reach a given accuracy in orbital elements. The prediction of future positions should be a stable procedure even in the presence of nearly-degenerate orbital elements. This will require that we fully understand the nature of the degeneracies from short arcs. We would like to find a parameterization of the orbits in which the degeneracies are confined to as few parameters as possible.
* The release of the USNO-A2.0 astrometric catalog111 The USNO-A2.0 catalog, produced by D. Monet et al., is described at http://www.nofs.navy.mil/projects/pmm/USNOSA2doc.html. now makes it possible to measure astrometric positions to 0$`\stackrel{}{\mathrm{.}}`$2 accuracy over the full sky (Deutsch, 1999) in a typical CCD image, and to produce reasonable estimates of the uncertainties on each position measurement. Uncertainty estimation can therefore proceed by the straightforward propagation of measurement errors. Space telescopes and adaptive-optics ground-based telescopes will commonly produce relative positions to 0$`\stackrel{}{\mathrm{.}}`$01 accuracy in the near future. Very accurate positions should in principle produce accurate orbits even over relatively short arcs.
* KBOs are at distances $`d30`$ AU so their apparent motions are dominated by reflex motion. The observed arcs are a small fraction of the orbital period $`P`$ even after a decade. Apparent motion patterns are thus very simple and unambiguous in comparison, say, to Near-Earth Objects.
* If we parameterize the distance to the object by $`\gamma =(1\mathrm{AU})/d0.03`$, the acceleration of the KBO is $`\gamma ^210^3`$ times smaller than Earth’s acceleration, and the component transverse to the line of sight is $`\gamma ^310^{4.5}`$ times Earth’s. The KBO motion is nearly inertial, and gravitational acceleration can be treated as a perturbation. Instead of parameterizing orbits by the usual element vector $`𝐚=\{a,e,i,\mathrm{\Omega },\omega ,T_p\}`$, an orbit is more stably specified by some Cartesian phase space vector $`𝐏=\{x_0,y_0,z_0,\dot{x}_0,\dot{y}_0,\dot{z}_0\}`$ at the time of initial observation.
* The annual parallax is limited to $`2\gamma 3\mathrm{°}`$, and the total object motion is only $`10\mathrm{°}`$ over a decade. We will therefore express sky positions as coordinates $`(\theta _x,\theta _y)`$ in a tangent-plane projection of the celestial sphere about an appropriate reference point, e.g. the first observed location. Furthermore, the line of sight to the object changes by only a few degrees, and we will find it convenient to align our Cartesian coordinates with the $`z`$ axis along the initial line of sight.
* Daily parallax is $`8\times 10^5\gamma 0\stackrel{}{\mathrm{.}}5`$. For most ground-based observations this will be too small to yield useful information.

### 2.2 Exact Equations

The input data are the observed right ascension and declination $`(\alpha _i,\delta _i)`$ of the target at times $`t_i`$. We will assume that all positions have been measured in the J2000 celestial reference frame via astrometry tied to background stars in the USNO-A2 catalog. Measuring relative to the background stars exactly cancels the effects of stellar aberration, and reduces the maximal error from gravitational deflection of starlight to the milli-arcsecond level, so we can treat the problem with Euclidean geometry. We transform the observed coordinates to a tangent-plane projection $`(\theta _x,\theta _y)`$ about a reference location $`(\alpha _0,\delta _0)`$. The $`\theta _x`$ axis is taken in the Eastward direction of the J2000 ecliptic coordinate system, and $`\theta _y`$ points toward J2000 ecliptic North. Formulae for the conversion from $`(\alpha ,\delta )`$ to $`(\theta _x,\theta _y)`$ are straightforward and summarized in Appendix A. Figure 1 illustrates the relation of our coordinate systems to the ecliptic.

We will generally take this projection axis to be the first observed position, and will take the first observation to be at time $`t_0=0`$. We will assume that $`\theta _x`$ and $`\theta _y`$ have equal and uncorrelated measurement uncertainties $`\sigma _i`$ (generalization to non-circular error regions is straightforward).

We define our inertial Cartesian coordinate system by placing the origin at (or near) the location of the observer at the initial epoch $`t=0`$. The $`z`$ axis is directed toward the reference direction $`(\alpha _0,\delta _0)`$ of the tangent-plane projection, i.e. along the initial line of sight. The $`y`$ axis is in the plane containing the $`z`$ axis and the J2000 North ecliptic pole, i.e. parallel to the $`\theta _y`$ axis of the tangent plane. The $`x`$ axis then lies in the ecliptic plane and points along $`\theta _x`$ to ecliptic East. Unless otherwise specified, we will take the units of distance and time to be AU and Julian year, respectively. The coordinate systems are illustrated in Figure 1.

If $`𝐱(t)`$ represents the position of the target at time $`t`$ and $`𝐱_E(t)`$ is the position of the observer on Earth, then the observed sky position at time $`t`$ is

$$\begin{array}{ccc}\theta _x(t)& =& \frac{x(t\mathrm{\Delta }t)x_E(t)}{z(t\mathrm{\Delta }t)z_E(t)}\\ \theta _y(t)& =& \frac{y(t\mathrm{\Delta }t)y_E(t)}{z(t\mathrm{\Delta }t)z_E(t)}\end{array}$$
(1)
where $`\mathrm{\Delta }t`$ is the light-travel time $`|𝐱(t\mathrm{\Delta }t)𝐱_E(t)|/c`$. The definition of $`\mathrm{\Delta }t`$ is recursive but in practice is easily solved for iteratively because of the lack of ambiguity in KBO orbit solutions.

We wish our uncertainties to be limited by the errors in the astrometric catalog, so our requirement is to model the target position to $`<0\stackrel{}{\mathrm{.}}1`$ accuracy over $`10`$-year time spans. For target distance $`30`$ AU, the object and observer positions must be correct to $`2000`$ km. The Earth’s position relative to the Sun and Solar System barycenter (SSB) is available from the JPL DE405 Ephemeris222 The DE405 ephemeris is available at http://ssd.jpl.nasa.gov/eph\_info.html (Standish et al., 1995) to an accuracy many orders of magnitude better than we require. A simplistic model of the Earth as an oblate spheroid in constant rotation gives the observatory position to accuracy well within our tolerance, so we may consider $`𝐱_E`$ to be known exactly. The interpolated JPL ephemeris and topocentric correction are rapidly calculable.

#### 2.2.1 Orbit Approximations

We write the target orbit as

$$𝐱(t)=𝐱_0+\dot{𝐱}_0t+𝐠(t),$$
(2)
where $`𝐠`$ is the gravitational perturbation satisfying

$`𝐠(t=0)`$ $`=`$ $`\dot{𝐠}(t=0)=0,`$ (3)
$`\ddot{𝐠}(t)`$ $`=`$ $`{\displaystyle \underset{j}{}}{\displaystyle \frac{GM_j[𝐱(t)𝐱_j(t)]}{|𝐱(t)𝐱_j(t)|^3}},`$ (4)
where the sum runs over the other bodies in the Solar System. We include forces from the Sun and the giant planets, obtaining their positions from the DE405 ephemeris. For objects beyond 30 AU, one could approximate the total gravitational field of the Solar System as a monopole originating the Solar System barycenter, with errors of only $`10^4`$ AU over 10 years (barring Neptune encounters), yielding positions accurate to $`<1\mathrm{}`$. The computational expense of the full N-body form Eq. (4) is insignificant, however, as the time steps in the orbit integration can be quite long (20 days or more) without significant loss of accuracy over a decade.

The combination of Eqs. (1), (2), and (4) define our model for the position $`(\theta _x,\theta _y)`$ as a function of time and the initial state vector $`𝐏`$. A comparison of our positions for Pluto, Neptune, and 1992 QB1 to those generated by the JPL ephemeris indicate residual differences of $`<0\stackrel{}{\mathrm{.}}05`$ for a 10-year integration.

### 2.3 Approximate Formulae & Degeneracies

We can understand the nature of the short-arc degeneracies by expanding Eq. (1) in powers of time $`t`$ and distance parameter $`\gamma `$. The formulae are simpler if we redefine our phase space basis in terms of the vector $`𝜶=\{\alpha ,\beta ,\gamma ,\dot{\alpha },\dot{\beta },\dot{\gamma }\}`$, where

$$\begin{array}{ccc}\alpha x_0/z_0\hfill & \beta y_0/z_0\hfill & \gamma 1/z_0\hfill \\ \dot{\alpha }\dot{x}_0/z_0\hfill & \dot{\beta }\dot{y}_0/z_0\hfill & \dot{\gamma }\dot{z}_0/z_0.\hfill \end{array}$$
(5)
Please note that $`\dot{\alpha }`$, $`\dot{\beta }`$, and $`\dot{\gamma }`$ quantify the motion along our three axes but are not the time derivatives of $`\alpha `$, $`\beta `$, and $`\gamma `$. We have chosen our coordinate system so that $`\alpha ,\beta 0`$, and we expect $`\dot{\alpha },\dot{\beta },\dot{\gamma }2\pi \gamma ^{3/2}`$ for bound orbits. If the orbit is nearly circular, the line-of-sight velocity is down another factor of $`\gamma `$ so that $`\dot{\gamma }2\pi \gamma ^{5/2}`$. The $`x`$ component of the apparent motion is now

$`\theta _x`$ $`=`$ $`{\displaystyle \frac{\alpha +\dot{\alpha }t+\gamma g_x(t)\gamma x_E(t)}{1+\dot{\gamma }t+\gamma g_z(t)\gamma z_E(t)}}`$ (6)
$``$ $`\alpha +(\dot{\alpha }\gamma \dot{x}_E)t`$
$`+\gamma (\ddot{g}_x\ddot{x}_E)t^2/2(\dot{\gamma }\gamma \dot{z}_E)(\dot{\alpha }\gamma \dot{x}_E)t^2`$
$`+O(t^3).`$
In the second line we have made the approximation that $`t1`$ year. The equation for $`\theta _y`$ is analogous. We see that the solution for the orbital parameters $`𝜶`$ has three regimes:

#### 2.3.1 Slope Regime

For arcs spanning consecutive nights, we can only determine the instantaneous positions $`\alpha `$ and $`\beta `$, plus the slopes

$`\dot{\theta }_x`$ $`=`$ $`(\dot{\alpha }\gamma \dot{x}_E)\dot{\alpha }\gamma (2\pi \mathrm{cos}\varphi _0)`$ (8)
$`\dot{\theta }_y`$ $`=`$ $`(\dot{\beta }\gamma \dot{y}_E)\dot{\beta }`$ (9)
The approximations assume that we are observing near the ecliptic plane, and that the Earth’s orbit is circular with phase $`\varphi _0`$ at time $`t=0`$ relative to the Sun-target vector. In this situation the line-of-sight motion $`\dot{\gamma }`$ is completely unconstrained, and there is a total degeneracy between the distance $`\gamma `$ and the ecliptic motion $`\dot{\alpha }`$. For a circular orbit, $`|\dot{\alpha }|(\sqrt{\gamma }/\mathrm{cos}\varphi _0)(\gamma \dot{x}_E)`$, and for any bound orbit $`|\dot{\alpha }|<(\sqrt{2\gamma }/\mathrm{cos}\varphi _0)(\gamma \dot{x}_E)`$. The usual observational strategy, therefore, is to observe near opposition ($`\varphi _0=0`$) to minimize the proper-motion term $`\dot{\alpha }`$ relative to the parallax term $`\gamma \dot{x}_E`$, giving an estimate of $`\gamma `$ to fractional accuracy of $`\sqrt{2\gamma }0.25`$ or better. An ephemeris may be predicted under the assumption of a prograde circular orbit. Below we will describe a simple way to place an error ellipse on the position without assuming a circular orbit.

#### 2.3.2 Acceleration Regime

As the arc length grows, the next bit of information gleaned is the apparent acceleration vector:

$`\ddot{\theta }_x`$ $`=`$ $`\gamma (\ddot{g}_x\ddot{x}_E)2(\dot{\gamma }\gamma \dot{z}_E)(\dot{\alpha }\gamma \dot{x}_E)`$ (10)
$``$ $`\gamma (\ddot{g}_x4\pi ^2\mathrm{sin}\varphi _0)+2\dot{\gamma }(2\pi \gamma \mathrm{cos}\varphi _0).`$ (11)
The second line again assumes an ecliptic observation and a circular Earth orbit, and keeps leading-order terms only. We see first that the effect of $`\dot{\gamma }`$ on the acceleration is a factor $`\gamma /2\pi <0.005`$ smaller than the $`\gamma \ddot{x}_E`$ term. In this regime, therefore, $`\dot{\gamma }`$ is still indeterminate. The target’s transverse acceleration $`\ddot{g}_x`$ is smaller by a factor $`\gamma ^3`$ than the Earth’s acceleration and can also be ignored at this point. The apparent acceleration is then

$`\ddot{\theta }_x`$ $``$ $`\gamma (4\pi ^2\mathrm{sin}\varphi _0)`$ (12)
$`\ddot{\theta }_y`$ $``$ $`0`$
We see that the apparent acceleration is a robust and simple way to constrain $`\gamma `$, the distance to the target. This method fails however, for observations at opposition ($`\varphi _0=0`$), where the reflex acceleration vanishes.

What time span is needed to measure the acceleration to a useful accuracy? Let us presume that the object is initially detected in a pair of observations separated by 24 hours, then recovered in a single observation some time $`T`$ later, and that each observation has astrometric uncertainty $`\sigma `$. Using Eq. (12) to determine $`\gamma `$, the fractional error in the distance to the target is then roughly

$`{\displaystyle \frac{\sigma _d}{d}}`$ $`=`$ $`{\displaystyle \frac{\sigma _\gamma }{\gamma }}{\displaystyle \frac{\sigma }{\sqrt{2}\pi ^2\gamma \mathrm{sin}\varphi _0(24\mathrm{h})T}}`$
$`=`$ $`7.5\%\left({\displaystyle \frac{\sigma }{0\stackrel{}{\mathrm{.}}2}}\right)\left({\displaystyle \frac{1\mathrm{week}}{T}}\right)\left({\displaystyle \frac{\mathrm{sin}45\mathrm{°}}{\mathrm{sin}\varphi _0}}\right)\left({\displaystyle \frac{d}{40\mathrm{AU}}}\right)`$
The distance to most KBOs can be determined from the apparent acceleration with only 1 week’s arc, as long as the observations are not near opposition. For a main belt asteroid, the apparent acceleration is an order of magnitude larger and is easily detected in a 24-hour span. This means it is easy to distinguish main belt asteroids from KBOs away from opposition, even near the main-belt turnaround points.

A typical (non-opposition) 1-week arc thus gives a useful estimate of 5 of the 6 phase-space parameters of the orbit, with the line-of-sight motion being completely indeterminate. In this basis the degeneracy is confined to a single parameter, whereas the degeneracy is shared between most of the traditional orbital elements. In particular neither $`a`$ nor $`e`$ is at all well constrained at this point.

It is interesting to note that the orbital angular momentum is sensibly constrained in this regime, since the line-of-sight velocity is nearly radial to the Sun and contributes little to the angular momentum. The energy of the orbit, however, is

$$E=GM_{}\gamma (12\gamma \mathrm{cos}\beta _0+\gamma ^2)^{1/2}+(\dot{\alpha }^2+\dot{\beta }^2+\dot{\gamma }^2)/2\gamma ^2,$$
(14)
and is uncertain until $`\dot{\gamma }`$ can be constrained. Here $`\beta _0`$ is the solar elongation of the target at the initial epoch.

#### 2.3.3 Fully Constrained Regime

The line-of-sight velocity can be constrained as the Earth moves around its orbit to give us a new point of view. From Eq. (6), the leading $`\dot{\gamma }`$ term in the sky position is $`\dot{\gamma }t\gamma x_E(t)\dot{\gamma }\gamma t`$. We know that $`\dot{\gamma }^2<2GM_{}\gamma ^3`$ from the simple constraint that the object be bound. A useful constraint on $`\dot{\gamma }`$ is obtained only after sufficient time that a $`\dot{\gamma }`$ value below the binding limit would produce a measured displacement that is distinguishable from terms in the other 5 parameters. This can occur on a timescale of months for high-quality observations. The uncertainty in $`\dot{\gamma }`$ (or $`\dot{z}`$) dominates the phase space uncertainties for essentially all lengths of arc, though it does not always dominate the ephemeris uncertainties.

### 2.4 Fitting Procedures

With 3 or more observations spanning several months or more, we will be in the fully-constrained regime. At this point the constraint of the orbital parameter vector $`𝜶`$ from the measured positions $`𝜽_i`$ becomes a straightforward non-linear minimization problem. The formalism for Bayesian orbit fitting is discussed in detail by Muinonen and Bowell (1993). We will assume Gaussian measurement errors, in which case the most likely parameter vector $`\widehat{𝜶}`$ is the least-squares value, minimizing

$$\chi ^2=\underset{i}{}\left\{\frac{[\theta _{x,i}\theta _x(\widehat{𝜶},t_i)]^2}{\sigma _i^2}+\frac{[\theta _{y,i}\theta _y(\widehat{𝜶},t_i)]^2}{\sigma _i^2}\right\}.$$
(15)
The model $`\theta _x(\widehat{𝜶},t_i)`$ for the predicted position is embodied by Eqs. (5), (1), (2), and (4). The minimization of $`\chi ^2`$ can use standard algorithms. We have implemented a slight modification of the Levenberg-Marquardt method from Press et al. (1988). For efficient location of the minimum, we require a good starting point for the non-linear search, and estimates of the derivative matrix $`d\theta /d𝜶`$.

The determination of the derivatives is again greatly simplified for distant objects. The gravitational acceleration is, over a $`10`$-year period, a small perturbation to the inertial motion, barring close encounters. The inertial motion will be $`2\pi \sqrt{\gamma }t`$ and the acceleration term will integrate to $`4\pi ^2\gamma ^2t^2/2`$. Furthermore the acceleration is toward the barycenter, so the transverse component is down by another factor of $`\gamma `$. The angular displacement due to gravity is only a few arcseconds per year, $`300`$ times smaller than the inertial motion and reflex motion. It is clear that in calculating the derivative matrix $`𝐱/𝜶`$ from Eq. (2) that we can ignore the gravitational term. All of the derivatives are thus quite trivial and there is no need to integrate the derivatives of the orbit.

Selection of a sensible starting point for the minimization is also simple for distant objects. We produce a version of the equations of motion which ignores gravity, is linearized, and also ignores the $`\dot{\gamma }`$ term:

$`\theta _x\alpha +\dot{\alpha }t\gamma x_E(t),`$ (16)
$`\theta _y\beta +\dot{\beta }t\gamma y_E(t),`$
These equations are linear in the parameter vector $`𝜶`$ and hence have a closed-form minimization requiring only the inversion of a $`5\times 5`$ matrix (Press et al., 1988). These very simple equations are actually good to a few arcseconds for the first year of observation, and in any case provide (with $`\dot{\gamma }=0`$) a good starting point for the non-linear minimization.

Using the analytic derivatives, the Levenberg-Marquardt algorithm also calculates the covariance matrix $`𝚺_\alpha `$ of the fitted parameters. Since the model $`𝜽(𝜶,t)`$ is dominated by terms linear in the parameters, the linearized approximation

$$\mathrm{\Delta }\chi ^2=(𝜶\widehat{𝜶})^T𝚺_\alpha (𝜶\widehat{𝜶})$$
(17)
will be accurate even for poorly constrained orbits. This would not have been the case had we chosen the orbital-element vector $`𝐚`$ as the parameter set.

### 2.5 Position and Orbital Element Estimation

Once the best-fit $`\widehat{𝜶}`$ vector and its covariance matrix are determined from the fit, determination of orbital elements and future positions is also simple. The predicted position at any time $`t`$ is $`𝜽(\widehat{𝜶},t)`$. The uncertainty ellipse for the position is specified by the covariance matrix

$$𝚺_\theta =\left(\frac{d𝜽(𝜶,t)}{d𝜶}\right)^T𝚺_\alpha \left(\frac{d𝜽(𝜶,t)}{d𝜶}\right).$$
(18)
The derivatives, as above, are rapidly calculated by ignoring gravity. Again the position is nearly linear in our chosen orbital basis, so this linearized transformation will be accurate even in the case of large uncertainties in the parameters.

The best-fit orbital element vector $`\widehat{𝐚}`$ is determined from the best-fit phase-space parameter set $`\widehat{𝜶}`$ via the usual equations, which are reviewed in Appendix B. The determination of an uncertainty ellipsoid in orbital-element space requires mapping of the covariance matrix with the derivatives $`d𝐚/d𝜶`$. Some of these derivatives are quite complicated analytic expressions, but can still be coded and rapidly evaluated. The map from $`𝜶`$ to $`𝐚`$ is non-linear, and in some cases degenerate, so the linearized covariance matrix is formally correct only when the uncertainties in the elements become small. We do, however, obtain a useful estimate of the element uncertainties even for less exact determinations.

#### 2.5.1 Fitting Singly Degenerate Orbits

If the arc is too short, the minimization of $`\chi ^2`$ will proceed but the best-fit parameters may be non-physical. The degeneracy in the “acceleration regime” will be manifested as one or more large elements of the covariance matrix. In particular, several-week orbits will have large values of $`\mathrm{\Sigma }_{\dot{\gamma }\dot{\gamma }}`$. We recognize these orbits as unrealistic because they may be unbound. The assumption of a bound orbit limits $`\dot{\gamma }`$ to (cf. Eq. )

$$\dot{\gamma }^2\dot{\gamma }_{\mathrm{bind}}^2=GM_{}\gamma ^3(1+\gamma ^22\gamma \mathrm{cos}\beta _0)^{1/2}\dot{\alpha }^2\dot{\beta }^2.$$
(19)
A conservative approach to assigning uncertainties in this case is to assume that $`\dot{\gamma }`$ has a uniform a priori probability in the interval $`[\dot{\gamma }_{\mathrm{bind}},+\dot{\gamma }_{\mathrm{bind}}]`$. If, however, we assign a Gaussian prior distribution to $`\dot{\gamma }`$, we can preserve our least-squares approach to the fitting problem. We therefore fit by assuming that $`\dot{\gamma }=0`$ with an RMS uncertainty of $`\sigma _{\mathrm{bind}}=\dot{\gamma }_{\mathrm{bind}}/\sqrt{3}`$. This value of $`\sigma _{\mathrm{bind}}`$ has the same variance as the uniform distribution, and the entire $`\pm \dot{\gamma }_{\mathrm{bind}}`$ range (plus some marginally unbound orbits) is contained within the $`2\sigma `$ contour, so we will obtain a realistic assessment of the uncertainties.

Our procedure for handling the $`\dot{\gamma }`$ degeneracy is thus:

1. Fit the orbit with all six parameters free to vary.
2. If the unconstrained fit yields $`\mathrm{\Sigma }_{\dot{\gamma }\dot{\gamma }}<\sigma _{\mathrm{bind}}^2`$, then $`\dot{\gamma }`$ is well constrained by the observations and we are done. Otherwise $`\dot{\gamma }`$ is essentially decoupled from the observations, so we proceed using the binding constraint instead.
3. Re-fit the observations but set $`\dot{\gamma }=0`$ and derive values and a covariance matrix for the remaining 5 parameters. Note that this best-fit orbit is very similar to the Väisälä orbit, as a line-of-sight velocity near zero implies an orbit near its perihelion or aphelion.
4. Augment the covariance matrix by assigning $`\mathrm{\Sigma }_{\dot{\gamma }\dot{\gamma }}=\sigma _{\mathrm{bind}}^2`$, placing zeros in the remaining off-diagonal elements.

The best-fit parameters $`\widehat{𝜶}`$ and covariance $`𝚺_\alpha `$ may then be used to predict positions and orbital elements as in the non-degenerate case. The $`2\sigma `$ error ellipses will encompass $`>95\%`$ of all bound orbits that are consistent with the data, and not be restricted to circular or perihelic orbits. If we have some preconceptions about the maximal ellipticity of the orbits, this may be incorporated into the position predictions by placing a tighter a priori constraint on $`\dot{\gamma }`$.

#### 2.5.2 Fitting Doubly Degenerate Orbits

With only a single or few nights’ data, or if there are only 2 observations, we are in the “slope regime.” The above technique will fail due the additional degeneracy between $`\dot{\alpha }`$ and $`\gamma `$ (and $`\dot{\beta }`$ for observations out of the ecliptic plane). This will manifest itself in the covariance matrix obtained in Step 3 above, the 5-parameter fit to the orbit. The uncertainties $`\mathrm{\Sigma }_{\dot{\alpha }\dot{\alpha }}`$ and/or $`\mathrm{\Sigma }_{\dot{\beta }\dot{\beta }}`$ will be large.

In this case we again introduce the binding constraint to limit the extent of the degeneracy. We express the transverse kinetic energy in terms of a parameter $`f_b`$:

$$\dot{\alpha }^2+\dot{\beta }^2=(1+f_b)GM_{}\gamma ^3.$$
(20)
The binding constraint is then $`1<f_b<1`$, with $`f_b=0`$ corresponding to a nearly circular orbit. We crudely implement this constraint by minimizing a new quantity

$$\stackrel{~}{\chi }^2=\chi ^2+f_b^2/3.$$
(21)
(with $`\chi ^2`$ from Eq. 15), which pushes the solution toward a circular orbit, while yielding a covariance matrix that reflects both the observational uncertainties and the possibility of non-circular orbits. More sophisticated constrained optimizations are possible but there is little point to seeking precise constraints on such a poorly determined orbit.

## 3 Implementation and Examples

### 3.1 C Code

All of the above algorithms have been implemented in a family of subroutines and drivers written in C. These include routines for: creation and access of binary versions of the DE405 ephemeris files (provided by David Hoffman); orbit integration; transformations among the various Cartesian and spherical coordinate systems; transformations between phase-space and orbital-element bases for orbits; reduction of positions to topocentric coordinates; and the fitting procedures described in the previous section. The code is portable and self-contained, and is available from the first author upon request.

The main program is fit\_radec, which takes as input a list of astrometric observations, and executes the orbit-fitting, including proper handling of degeneracies. The output is a file specifying the best-fit $`𝜶`$ description of the orbit, and its full covariance matrix $`𝚺_\alpha `$.

A program abg\_to\_aei transforms the $`𝜶`$ parameters into traditional orbital elements. The covariance matrix is also transformed with a linear approximation.

The $`𝜶`$ information may be passed to the program predict, which generates predicted positions and their error ellipses for any specified observation date and site.

The entire process is extremely fast on the typical present-day desktop computer. A fit to 275 mock observations of Pluto over a 14-year arc takes about 1 second. The entire ensemble of observations of all $`300`$ KBOs reported to the MPC are fit in a few seconds.

### 3.2 Example and Verification

We have verified the performance of the algorithms and software by conducting mock observations of Pluto, Neptune, and several KBOs. Figure 2 shows the results of an ensemble of simulated observations of Pluto from 1995–1998. The true position of Pluto is obtained from the DE405 ephemeris; mock observations are created at specified dates by adding Gaussian deviations to the true positions. An orbit is then fit to the mock observations. This fitted orbit is used to estimate the osculating barycentric elements $`a`$ and $`e`$ for Pluto as well as the uncertainty ellipse in the $`a`$$`e`$ plane. Alternatively we may use the orbit fit to estimate the $`(\theta _x,\theta _y)`$ sky position of Pluto at some future date and its projected error ellipse. We can then examine whether the estimated orbital elements or future position are consistent with the true ephemeris values. In all cases we find the mean $`\chi ^2`$ value of the fit to the orbit to be as expected from the degrees of freedom in the fit.

The first panel in Figure 2 shows the results of using 4 observations over an 11-day arc (JD 2449927–2449938) to predict the location of Pluto $`300`$ days after the first observation. Each point on the plot is the prediction from one of 1000 realizations of the mock 11-day arc. The star shows the true location of Pluto at the 300-day time. The ellipses show the $`1\sigma `$ and $`2\sigma `$ contours of the estimated error on the position given the observational errors. The estimated uncertainty ellipse is quite consistent with the scatter in the mock predictions. There appear to be two flaws in our prediction: first, the mean of the predicted positions does not coincide with the true position; second, the minor axis of the uncertainty ellipse appears to be significantly larger than the scatter. Both of these phenomena are expected, however, because the 11-day arc is the “acceleration regime” and is degenerate in $`\dot{\gamma }`$. The fits therefore have assumed a value and uncertainty for $`\dot{\gamma }`$. The displacement of the simulations from the true value is due to our assumption that $`\dot{\gamma }=0`$ for each fit, which is not the true value. But the uncertainty ellipse includes a contribution from our ignorance of $`\dot{\gamma }`$, which makes the ellipse extend beyond the measured scatter, and in fact properly include the true position.

In the second panel we simulate an observation sequence consisting of the 11-day arc plus and additional observation at 300 days. We plot the estimated $`a`$ and $`e`$ values from 1000 mock observations of the 300-day arc. In this case the uncertainty ellipses give a roughly correct estimate of the size of the errors, but the scatter does not match the error ellipses in detail. This demonstrates a limitation of our method, which is that the propagation of the covariance matrix $`𝚺_\alpha `$ to the orbital elements is incorrect if the map between them is significantly nonlinear. This is the case here, as Pluto is very near perihelion at this time. The 300-day arc determines the distance to Pluto quite well, and hence the scatter plot is bounded by a lower limit to the perihelion $`a(1e)`$.

The third panel shows the results of using the 300-day arc to predict the position after 1000 days, and the fourth panel shows the scatter in the $`a`$$`e`$ plane of elements determined from 6 observations on a 1000-day arc. We can now usefully constrain all 6 parameters of the orbit, and the scatter about the true position and orbit are very accurately described by the estimated covariance matrices.

## 4 Optimal Strategies

An optimal schedule for recovery of KBOs would maximize some measure of the accuracy of the KBO orbit while minimizing the number of recovery observations required. A constraint on the observing schedule is that the positional uncertainty must not become larger than the field of view of the instrument, lest the recovery be missed and the object be lost. As an illustration of the utility of the methods described in this paper, we derive an optimal schedule of observations using a crude algorithm.

We take our measure of the quality of the orbit fit to be the uncertainty $`\sigma _a=\sqrt{\mathrm{\Sigma }_{aa}}`$ in the semi-major axis of the orbit. This is not meant to be a universal metric for quality but rather a simple example. We place two constraints on the observing schedule: first, observations are not allowed when the solar elongation of the target is less than 90°. Second, we take the field of view of the recovery instrument to be $`10\mathrm{}`$; then the 1-sigma positional uncertainty must be $`2\stackrel{}{\mathrm{.}}5`$ if the full 2-sigma error ellipse is to be contained within the field of view for $`>95\%`$ chance of recovery. Note that many CCD mosaics currently on 2-meter and 4-meter class telescopes have fields of view that substantially exceed 10′.

A rigourous optimization for a set of $`N`$ observations would consider all possible sets of $`N`$ nights and choose the set which minimizes the final $`\sigma _a`$. We adopt a somewhat simpler “greedy” algorithm for optimization as follows: after observation $`M`$, we consider all nights for which observation $`M+1`$ would be possible, namely those satisfying the elongation ad positional uncertainty costraints. We next calculate the $`\sigma _a`$ which would result from an observation on each possible night; the $`M+1`$ observation is then scheduled for the night which would result in minimal $`\sigma _a`$.

Figure 3 illustrates the observation schedule derived for hypothetical discovery of Pluto in the year 1995. We assume that the object is discovered on a pair of nights roughly 45 days before opposition, then confirmed on another pair of nights 10 days later. We also presume that each observation has an error of $`\sigma =0\stackrel{}{\mathrm{.}}2`$. As seen in the Figure, the object must be observed once again in the discovery season as its uncertainty is too high to be recovered when it emerges from behind the Sun the following year. The object is now sufficiently well characterized to be found for nearly 3 years. When is the best time during this window to observe again? Interestingly the accuracy on $`a`$ improves only weakly the longer we wait for our next observation. The “greedy” strategy suggests we wait 3 years, when the positional uncertainty approaches 2$`\stackrel{}{\mathrm{.}}`$5. So after the initial 4 discovery/confirmation images and only 2 recoveries, $`\sigma _a0.03`$ AU, and the orbit is known well enough to permit recovery for the next decade.

A counterintuitive result is that the timing of further observations has little effect upon the accuracy of $`a`$—only increased number and/or accuracy of observations matters at this point. This is because the uncertainty in $`a`$ derives primarily from the uncertainty in $`\gamma `$, i.e. the distance to the target. The leading observable consequence of $`\gamma `$ is the parallax $`\gamma x_E(t)`$, which does not grow with time beyond 1 year—hence the very weak improvement with longer arcs. Indeed we find that observing near the two quadratures of the second year yields a $`\sigma _a`$ that is as low as spacing these two observations over several years.

Another common observing scheme is to run discovery observations during two nights near opposition, then attempt confirmation one month later. For Pluto, the positional error is sufficiently small this second month to permit reliable confirmation. Henceforth the optimal strategy plays out just as the previous scenario: a further observation is required before quadrature to keep the object from being lost while behind the Sun, and the orbital quality depends primarily upon the number, not the arc length, of recoveries in later years.

If the astrometric accuracy of the initial 11-day arc can be improved to 0$`\stackrel{}{\mathrm{.}}`$1, and the recovery tolerance raised to $`\sigma _a10\mathrm{}`$, then no further observations are required in the first year. The object is recoverable at the beginning of the second season, two observations (near the quadratures) in this second year reduce $`\sigma _a`$ to 0.0002 AU, and the position is good for a decade or more.

These results are of course dependent upon the orbit of the object and upon the nature of one’s metric for orbital quality. Our optimization method is, however, easily adapted to other metrics.

## 5 Summary

Orbit fits for KBOs and other distant objects are greatly simplified by the fact that the gravitational perturbation to the orbit is small compared both to the inertial motion and to the acceleration of the Earth. Gravity can be treated as a perturbation even for decades-long arcs, which simplifies the fitting process and the propagation of errors because there is no need to integrate the derivatives of the orbit. The orbit itself can be integrated quickly because large time steps can be taken. Distant objects move slowly across the sky, so we can also use the convenience of a tangent-plane projection of the celestial sphere. The dominance of reflex motion over proper motion means that KBO orbit solutions are unambiguous.

Exploiting these properties, we have created algorithms and software which quickly and accurately calculates orbital elements and ephemerides and their associated uncertainties for targets $`10`$ AU from the Sun. In cases where the observed arc is short enough to leave the orbit degenerate, we can still calculate sensible positional errors because we have chosen an orbital basis in which the degeneracy is confined to one parameter, namely the line-of-sight velocity. For very short arcs in the “slope regime,” 2 additional parameters become degenerate. We are still able to place error ellipses on ephemerides which are conservative in the sense that they will include any bound orbit. Most short-arc objects will be recovered closer to the MPC predicted ephemeris than our error ellipses would suggest, because the philosophy of the MPC is to use the stronger a priori consideration that the orbits resemble those of the majority of the known populations.

Our routines should prove valuable for planning of recovery observations, as we have illustrated with some simple examples. Estimates of the uncertainty ellipsoids in orbital-element space will also be needed for study of the long-term dynamics of the known KBOs.

In the near future we expect that KBO positions will be measured to accuracies of a few milli-arcseconds with orbiting observatories—though such observations will be even more expensive than present-day KBO observations. It will then be possible to place strong constraints on KBO orbits even from very short (e.g. 24-hour) arcs. For example, the motion of an observatory in low-Earth orbit imparts a reflex motion on the KBO that is easily detected at this level, giving a measure of the distance to the target in less than an hour. The methods described herein will provide the means to plan and to analyze such observations.

We thank David Hoffman for making available his C code to read the DE405 ephemerides; Brian Marsden for discussions on fitting of KBO orbits; Julia Plummer for assistance with Figure 1; and Renu Malhotra and Lynne Allen for serving as “guinea pigs” for tests of the software. This work is supported by NASA Planetary Astronomy grant #NAG5-7860, and grant #AST-9624592 from the National Science Foundation.

## Appendix A Coordinate Transformations

### A.1 Angular Coordinates

We make use of three sets of coordinates on the celestial sphere: equatorial coordinates $`(\alpha ,\delta )`$ in the ICRS reference frame, J2000; ecliptic latitude and longitude $`(\mathrm{},b)`$, also taken in the epoch 2000 frame; and our projected coordinates $`(\theta _x,\theta _y)`$ measured in a tangent plane about some reference direction $`(\alpha _0,\delta _0)`$ or $`(\mathrm{}_0,b_0)`$. The transformation from equatorial to ecliptic coordinates is implemented as

$`\mathrm{sin}b`$ $`=`$ $`\mathrm{cos}ϵ\mathrm{sin}\delta \mathrm{sin}ϵ\mathrm{cos}\delta \mathrm{sin}\alpha `$ (A1)
$`\mathrm{tan}\mathrm{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{cos}ϵ\mathrm{cos}\delta \mathrm{sin}\alpha +\mathrm{sin}ϵ\mathrm{sin}\delta }{\mathrm{cos}\delta \mathrm{cos}\alpha }},`$ (A2)
where $`ϵ=23\stackrel{}{\mathrm{.}}43928`$ is the obliquity of the ecliptic at J2000 (Standish et al., 1995). The inverse transformation simply reverses the sign of $`ϵ`$.

The map from ecliptic to tangent-plane coordinates about $`(\mathrm{}_0,b_0)`$ is

$`\theta _x`$ $`=`$ $`{\displaystyle \frac{\mathrm{cos}b\mathrm{sin}(\mathrm{}\mathrm{}_0)}{\mathrm{sin}b_0\mathrm{sin}b\mathrm{cos}b_0\mathrm{cos}b\mathrm{cos}(\mathrm{}\mathrm{}_0)}}`$ (A3)
$`\theta _y`$ $`=`$ $`{\displaystyle \frac{\mathrm{cos}b_0\mathrm{cos}b\mathrm{sin}b_0\mathrm{sin}b\mathrm{cos}(\mathrm{}\mathrm{}_0)}{\mathrm{sin}b_0\mathrm{sin}b+\mathrm{cos}b_0\mathrm{cos}b\mathrm{cos}(\mathrm{}\mathrm{}_0)}}.`$ (A4)
The inverse transformation is

$`\mathrm{sin}b`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}b_0+\theta _y\mathrm{cos}b_0}{\sqrt{1+\theta _x^2+\theta _y^2}}}`$ (A5)
$`\mathrm{sin}(\mathrm{}\mathrm{}_0)`$ $`=`$ $`{\displaystyle \frac{\theta _x/\mathrm{cos}b_0}{\sqrt{1+\theta _x^2+\theta _y^2}}}.`$ (A6)
Partial derivative matrices are straightforward for all transformations, and approximations for $`\theta _x,\theta _y1`$ are particularly simple.

### A.2 Spatial Coordinates

There are also three relevant spatial coordinate systems. The equatorial $`𝐱_{\mathrm{eq}}`$ has origin at the Solar System barycenter, with $`x`$ positive to the Vernal Equinox and $`z`$ positive to the North equatorial pole of J2000. The ecliptic system $`𝐱_{\mathrm{ec}}`$ also has origin at the barycenter, but $`z`$ is positive to the North ecliptic pole of J2000. Our “telescope-centric” coordinate system $`𝐱_T`$ has origin at the location of the observer at $`t=0`$, with $`z`$ positive along the line of sight to the target at $`t=0`$, the $`y`$ axis in the plane of $`z`$ and the North ecliptic pole—see Figure 1.

The transformation from equatorial to ecliptic coordinates is

$$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]_{\mathrm{ec}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}ϵ& \mathrm{sin}ϵ\\ 0& \mathrm{sin}ϵ& \mathrm{cos}ϵ\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]_{\mathrm{eq}}$$
(A7)
The inverse transformation inverts the sign of $`ϵ`$.

The transformation from ecliptic to telescope-centric coordinates is

$$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]_T=\left[\begin{array}{ccc}\mathrm{sin}\mathrm{}_0& \mathrm{cos}\mathrm{}_0& 0\\ \mathrm{cos}\mathrm{}_0\mathrm{sin}b_0& \mathrm{sin}\mathrm{}_0\mathrm{sin}b_0& \mathrm{cos}\mathrm{}_0\\ \mathrm{cos}\mathrm{}_0\mathrm{cos}b_0& \mathrm{sin}\mathrm{}_0\mathrm{cos}b_0& \mathrm{sin}\mathrm{}_0\end{array}\right]\left[\begin{array}{c}xx_0\\ yy_0\\ zz_0\end{array}\right]_{\mathrm{ec}}.$$
(A8)
The ecliptic-coordinate location of the observer at time $`t=0`$ is given by the vector $`\{x_0,y_0,z_0\}`$. The inverse transformation is

$$\left[\begin{array}{c}xx_0\\ yy_0\\ zz_0\end{array}\right]_{\mathrm{ec}}=\left[\begin{array}{ccc}\mathrm{sin}\mathrm{}_0& \mathrm{cos}\mathrm{}_0\mathrm{sin}b_0& \mathrm{cos}\mathrm{}_0\mathrm{cos}b_0\\ \mathrm{cos}\mathrm{}_0& \mathrm{sin}\mathrm{}_0\mathrm{sin}b_0& \mathrm{sin}\mathrm{}_0\mathrm{cos}b_0\\ 0& \mathrm{cos}b_0& \mathrm{sin}b_0\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]_T.$$
(A9)
Since all the transformations are linear, the partial derivatives are trivial.

## Appendix B Orbit Basis Transformations

The orbit-fitting algorithms produce a description of the orbit in terms of $`𝜶=\{\alpha ,\beta ,\gamma ,\dot{\alpha },\dot{\beta },\dot{\gamma }\}`$. To obtain orbital elements, we first transform to Cartesian phase space elements in the telescope-centric coordinate system:

$$\begin{array}{ccc}x_T=\alpha /\gamma \hfill & y_T=\beta /\gamma \hfill & z_T=1/\gamma \hfill \\ \dot{x}_T=\dot{\alpha }/\gamma \hfill & \dot{y}_T=\dot{\beta }/\gamma \hfill & \dot{z}_T=\dot{\gamma }/\gamma \hfill \end{array}$$
(B1)
The telescope-centric phase space coordinates can then be transformed to barycentric ecliptic Cartesian phase-space coordiates using Eq. (A9). The transformation matrix for the velocities is the same as that for the position, but there is no telescope-to-barycenter offset for the velocities.

With the barycentric position and velocity $`𝐱`$ and $`𝐯`$ in hand, the osculating orbital elements are determined with the standard formulae below. Using the gravitational constant $`\mu GM`$, we have

Semimajor axis:

$$a^1=\frac{2}{x}\frac{v^2}{\mu }.$$
(B2)
Eccentricity:

$`𝐞`$ $``$ $`\left[{\displaystyle \frac{v^2}{\mu }}{\displaystyle \frac{1}{x}}\right]𝐱\left[{\displaystyle \frac{𝐱𝐯}{\mu }}\right]𝐯`$ (B3)
$`e`$ $`=`$ $`|𝐞|`$ (B4)
With the further definitions

$`𝐡`$ $``$ $`𝐱\times 𝐯`$ (B5)
$`𝐧`$ $``$ $`\widehat{𝐳}\times 𝐡`$ (B6)
the inclination, ascending node, and argument of perihelion are given by

$`\mathrm{cos}i`$ $`=`$ $`{\displaystyle \frac{h_z}{h}}`$ (B7)
$`\mathrm{cos}\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \frac{n_x}{n}}`$ (B8)
$`\mathrm{cos}\omega `$ $`=`$ $`{\displaystyle \frac{𝐧𝐞}{ne}}.`$ (B9)
The eccentric and mean anomalies are

$`\mathrm{tan}E`$ $`=`$ $`{\displaystyle \frac{\overline{y}/b}{(\overline{x}/a)+e}}`$ (B10)
$`M`$ $`=`$ $`Ee\mathrm{sin}E,`$ (B11)
with

$`\overline{x}`$ $`=`$ $`{\displaystyle \frac{px}{e}}`$ (B12)
$`\overline{y}`$ $`=`$ $`{\displaystyle \frac{𝐱𝐯}{e}}\sqrt{{\displaystyle \frac{p}{\mu }}}`$ (B13)
$`b`$ $`=`$ $`a\sqrt{1e^2}`$ (B14)
$`p`$ $`=`$ $`{\displaystyle \frac{h^2}{\mu }}`$ (B15)
Finally, the time of periapse passage is

$$T_p=t_0M\sqrt{a^3/\mu },$$
(B16)
where $`t_0`$ is the time at which $`𝐱`$ and $`𝐯`$ are determined.

The partial derivatives necessary for propagation of errors are tedious but calculable. Our phase-space elements are centered on the solar system barycenter, therefore the osculating elements are also barycentric. We take $`M`$ as total mass of the solar system since we are using the barycenter.


text:
# World-Wide Web scaling exponent from Simon’s 1955 model

\[

## Abstract

Recently, statistical properties of the World-Wide Web have attracted considerable attention when self-similar regimes have been observed in the scaling of its link structure. Here we recall a classical model for general scaling phenomena and argue that it offers an explanation for the World-Wide Web’s scaling exponent when combined with a recent measurement of internet growth.

today \]

A quantity important for searching the World-Wide Web  is the number $`k`$ of links that point to a particular web page. Its probability distribution $`P(k)`$ exhibits power-law scaling  $`P(k)k^\gamma `$ that is not readily explained by standard random graph theory . An elegant model for scaling in copy and growth processes was proposed by Simon  in 1955 which describes scaling behaviour as observed in distributions of word frequencies in texts or population figures of cities . It models the dynamics of a system of elements with associated counters (e.g., words and their frequencies in texts, or nodes in a network and their connectivity $`k`$) where the dynamics of the system is based on constant growth via addition of new elements (new instances of words) as well as incrementing the counters (new occurrences of a word) at a rate proportional to their current values.

Reformulating this to model network growth consider a network with $`n`$ nodes with connectivities $`k_i`$, $`i=1\mathrm{}n`$, forming classes $`[k]`$ of $`f(k)`$ nodes with identical connectivity $`k`$. Iterate the following steps:

(i) With probability $`\alpha `$ add a new node and attach a link to it from an arbitrarily chosen node.

(ii) Else add one link from an arbitrary node to a node $`j`$ of class $`[k]`$ chosen with probability

$$P_{\text{new link to class }[k]}kf(k).$$
(1)
For this stochastic process, Simon finds a stationary solution exhibiting power-law scaling with exponent

$$\gamma =1+\frac{1}{1\alpha }.$$
(2)
The only free parameter of the model $`\alpha `$ reflects the relative growth of number of nodes versus number of links. In general small values of $`\alpha `$, therefore, predict scaling exponents near $`\gamma 2`$.

Let us apply this process to model the evolution of the World-Wide Web, identifying nodes with web pages. Data from two recent comprehensive Altavista crawls  provide an estimate for $`\alpha `$ in the present internet. These two measurements counted 203 million pages and 1466 million links in May 1999, and 271 million pages and 2130 million links in October 1999. The probability for adding a new web page is estimated from the observed increase in counts to $`\alpha 0.10`$. The subsequent prediction of Simon’s model for the exponent of the link distribution is $`\gamma =2.1`$ comparing well to current experimental results $`\gamma =2.1\pm 0.1`$  and $`\gamma =2.09`$ .

To compare with recently proposed models it may be interesting to note that the model by Barabási and Albert  can be mapped to the subclass $`\alpha =1/2`$ of Simon’s model, when using the simpler probability for a node being connected to another node $`i`$ with connectivity $`k_i`$

$$P_{\text{new link to i}}k_i.$$
(3)
Note that (3) implies (1) whereas the reverse is not true. Otherwise both models are based on the same two assumptions of growth and preferential linking. From this viewpoint, it is insightful to reconsider a recent discussion of their model. Adamic and Huberman point out that the “rich-get-richer” behaviour of single nodes imposed by (3) correlates age and connectivity of nodes . This, however, is disproven by data they present. They suggest to (and Barabási et al. in response show how to ) add individual growth rates to each node. While this solves the correlation problem, the price to pay is a large number of free parameters in the extended model. A simple solution to this problem has already been provided by Simon: Linking is guided by (1) instead of (3), considering not single nodes but classes of nodes with identical connectivities. This allows for different growth rates among class members, leaving just one free parameter. Above we determine this parameter from experimental data, enabling Simon’s classical scaling model to estimate the connectivity exponent of the World-Wide Web to $`\gamma =2.1`$.