text:
# RELATIVISTIC EPICYCLES: another approach to geodesic deviations

## 1 Introduction

The problem of motion of planets in General Relativity, considered as test particles moving along geodesic lines in the metric of Schwarzschild’s solution, has been solved in an approximate way by Einstein , who found that the perihelion advance during one revolution is given in the near-Keplerian limit by the formula

$$\mathrm{\Delta }\varphi =\frac{6\pi GM}{a(1e^2)}$$
(1)
where $`G`$ is Newton’s gravitational constant, $`M`$ the mass of the central body, $`a`$ the greater half-axis of planet’s orbit and $`e`$ its eccentricity.
This formula is deduced from the exact solution of the General Relativistic problem of motion of a test particle in the field of Schwarzschild metric, which leads to the expression of the angular variable $`\phi `$ as an elliptic integral, which is then evaluated after expansion of the integrand in terms of powers of the small quantity $`\frac{GM}{r}`$.
The formula has been successfully confronted with observation, giving excellent fits not only for the orbits with small eccentricities (e.g., one of the highest values of $`e`$ displayed by the orbit of Mercury, is $`e=0.2056`$ ), but also in the case when $`e`$ is very high, as for the asteroid Icarus ($`e=0.827`$), and represents one of the best confirmations of Einstein’s theory of gravitation. In the case of small eccentricities the formula (1) can be developed into a power series:

$$\mathrm{\Delta }\varphi =\frac{6\pi GM}{a}(1+e^2+e^4+e^6+\mathrm{}).$$
(2)
One can note at this point that even for the case of planet Mercury, the series truncated at the second term, i.e., taking into account only the factor $`(1+e^2)`$ will lead to the result that differs only by $`0.18\%`$ from the result predicted by relation (1), which is below the actual error bar.
This is why we think it is useful to present an alternative way of treating this problem, which is based on the use of geodesic deviation equations of first and higher orders. Instead of developing the exact formulae of motion in terms of powers of the parameter $`\frac{GM}{r}`$, we propose to start with an exact solution of a particularly simple form (i.e. a circular orbit with uniform angular velocity), and then generate the approximate solutions as geodesics being close to this orbit.
One of the advantages of this method is the fact that it amounts to treating consecutively systems of linear equations with constant coefficients, all of them being of harmonic oscillator type, eventually with an extra right-hand side being a known periodic function of the proper time. The approximate solution obtained in this manner has the form of a Fourier series and represents the closed orbit as a superposition of epicycles with diminishing amplitude as their circular frequencies grow as multiples of the basic one. This approach is particularly well-suited for using numerical computations. An example is provided by the computation of gravitational radiation from non-circular orbits, for which we use the well-known quadrupole formula.

## 2 Geodesic deviations of first and higher orders

Of many equivalent derivations of geodesic deviation equation we present the one which most directly leads to the results used in subsequent applications. Given a (pseudo)-Riemannian manifold $`V_4`$ with the line element defined by metric tensor $`g_{\mu \nu }(x^\lambda )`$,

$$ds^2=g_{\mu \nu }(x^\lambda )dx^\mu dx^\nu ,$$
(3)
a smooth curve $`x^\lambda (s)`$ parametrized with its own length parameter (or proper time) $`s`$ is a geodesic if its tangent vector $`u^\mu =(dx^\mu /ds)`$ satisfies the equation:

$$u^\lambda _\lambda u^\mu =0\frac{Du^\mu }{Ds}=\frac{du^\mu }{ds}+\mathrm{\Gamma }_{\lambda \rho }^\mu u^\lambda u^\rho =0.$$
(4)
where $`\mathrm{\Gamma }_{\rho \lambda }^\mu `$ denote the Christoffel connection coefficients of the metric $`g_{\mu \nu }`$.
Suppose that a smooth congruence of geodesics is given, of which the geodesics are labeled by a continuous parameter $`p`$: $`x^\mu =x^\mu (s,p),`$ such that the two independent tangent vector fields are defined by:

$$u^\mu (s,p)=\frac{x^\mu }{s}\mathrm{and}n^\mu (s,p)=\frac{x^\mu }{p}.$$
(5)
It is easily established that the rates of change of the tangent vectors in the mutually defined directions are equal:

$$n^\lambda _\lambda u^\mu =u^\lambda _\lambda n^\mu \frac{Du^\mu }{Dp}=\frac{Dn^\mu }{Ds}=\frac{^2x^\mu }{ps}+\mathrm{\Gamma }_{\lambda \rho }^\mu \frac{x^\lambda }{p}\frac{x^\rho }{s},$$
(6)
by virtue of the symmetry of Christoffel symbols in their lower indices.
The Riemann tensor can be defined using covariant derivations along the two independent directions of the congruence:

$$[u^\lambda _\lambda ,n^\rho _\rho ]Y^\mu =\left[\frac{D}{Ds}\frac{D}{Dp}\frac{D}{Dp}\frac{D}{Ds}\right]Y^\mu =R_{\lambda \rho \sigma }^\mu \frac{x^\lambda }{s}\frac{x^\rho }{p}Y^\sigma .$$
(7)
Replacing $`Y^\mu `$ by $`u^\mu `$ in the above formula, we get

$$[u^\lambda _\lambda ,n^\rho _\rho ]u^\mu =R_{\lambda \rho \sigma }^\mu \frac{x^\lambda }{s}\frac{x^\rho }{p}u^\sigma =R_{\lambda \rho \sigma }^\mu u^\lambda u^\sigma n^\rho .$$
(8)
By virtue of the geodesic equation (4) and Eq. (6), this can be written as

$$u^\lambda _\lambda \left(n^\rho _\rho u^\mu \right)=\frac{D}{Ds}\frac{Du^\mu }{Dp}=\frac{D^2n^\mu }{Ds^2}=R_{\lambda \rho \sigma }^\mu u^\lambda u^\sigma n^\rho .$$
(9)
This first-order geodesic deviation equation is often called the Jacobi equation, and is manifestly covariant.
In certain applications, Eq. (9) can be replaced by its more explicit, although non-manifestly covariant version:

$$\frac{d^2n^\mu }{ds^2}+2\mathrm{\Gamma }_{\lambda \rho }^\mu u^\lambda \frac{dn^\rho }{ds}+_\sigma \mathrm{\Gamma }_{\lambda \rho }^\mu u^\lambda u^\rho n^\sigma =0.$$
(10)
In this form of the geodesic deviation equation one easily identifies the relativistic generalizations of the Coriolis-type and centrifugal-type inertial forces, represented respectively by the second and third terms of Eq. (10).
The geodesic deviation can be used to construct geodesics $`x^\mu (s)`$ close to a given reference geodesic $`x_0^\mu (s)`$, by an iterative procedure as follows. Let the two geodesics be members of a congruence as above, with

$$x^\mu (s)=x^\mu (s,p),x_0^\mu (s)=x^\mu (s,p_0).$$
(11)
It follows by direct Taylor expansion, that

$`x^\mu (s,p)`$ $`=`$ $`x^\mu (s,p_0)+(pp_0){\displaystyle \frac{x^\mu }{p}}|_{(s,p_0)}+{\displaystyle \frac{1}{2}}(pp_0)^2{\displaystyle \frac{^2x^\mu }{p^2}}|_{(s,p_0)}+\mathrm{}`$ (12)
$`=`$ $`x_0^\mu (s)+\delta x^\mu (s)+{\displaystyle \frac{1}{2!}}\delta ^2x^\mu (s)+{\displaystyle \frac{1}{3!}}\delta ^3x^\mu (s)+\mathrm{},`$
where the more compact notation $`\delta ^nx^\mu (s)`$ describes the $`n^{th}`$-order geodesic deviation. Because $`(pp_0)`$ is supposed to be a small quantity, for convenience we may denote it $`ϵ`$. The first-order deviation is a vector, $`\delta x^\mu (s)=(pp_0)n_0^\mu (s)=ϵn_0^\mu (s)`$. But the second-order deviation is not a vector, and is given by

$$\delta ^2x^\mu (s)=(pp_0)^2\left(b^\mu \mathrm{\Gamma }_{\lambda \nu }^\mu n^\lambda n^\nu \right)_0=ϵ^2\left(b^\mu \mathrm{\Gamma }_{\lambda \nu }^\mu n^\lambda n^\nu \right)_0$$
(13)
where the covariant second-order deviation vector $`b^\mu `$ is defined by

$$b^\mu =\frac{Dn^\mu }{Dp}=\frac{n^\mu }{p}+\mathrm{\Gamma }_{\lambda \nu }^\mu n^\lambda n^\nu .$$
(14)
Straightforward covariant differentiation of Eq. (9), plus use of the Bianchi and Ricci identities for the Riemann tensor, implies that this second-order deviation vector $`b^\mu (s)`$ satisfies an inhomogeneous extension of the first-order geodesic deviation equation:

$$\frac{D^2b^\mu }{Ds^2}+R_{\rho \lambda \sigma }^\mu u^\lambda u^\sigma b^\rho =[_\nu R_{\lambda \rho \sigma }^\mu _\lambda R_{\nu \sigma \rho }^\mu ]u^\lambda u^\sigma n^\rho n^\nu +4R_{\lambda \rho \sigma }^\mu u^\lambda n^\rho \left(\frac{Dn^\sigma }{Ds}\right).$$
(15)
A more detailed formal derivation of this equation is given in appendix 2.

A rigorous mathematical study of geodesic deviations up to the second-order, as well as geometric interpretation, but using a different derivation, was presented in Ref. . Also, a Hamilton–Jacobi formalism has been derived in Refs. , which was applied to the problem of free falling particles in the Schwarzschild space-time . Fine effects resulting from the analysis of geodesic deviations of test particles suspended in hollow spherical satellites have been discussed in Ref. .

Obviously the procedure can be extended to arbitrarily high order geodesic deviations $`\delta ^nx^\mu (s)`$. This is of considerable practical importance, as it allows to construct a desired set of geodesics in the neighborhood of the reference $`x_0^\mu (s)`$, when the congruence of geodesics is not given a priori in closed form. Indeed, all that is needed is the set of deviation vectors $`(n_0^\mu (s),b_0^\mu (s),\mathrm{})`$ on the reference geodesic; obviously these vectors are completely specified as functions of $`s`$ by solving the geodesic deviation equations (9), (15) and their extensions to higher order, for given $`x_0^\mu (s)`$.
As in the case of the first-order deviation, it is sometimes convenient to write equation (15) in the equivalent but non-manifest covariant form

$$\begin{array}{c}\frac{d^2b^\mu }{ds^2}+_\rho \mathrm{\Gamma }_{\lambda \sigma }^\mu u^\lambda u^\sigma b^\rho +2\mathrm{\Gamma }_{\lambda \sigma }^\mu u^\lambda \frac{db^\sigma }{ds}=4\left(_\lambda \mathrm{\Gamma }_{\sigma \rho }^\mu +\mathrm{\Gamma }_{\sigma \rho }^\nu \mathrm{\Gamma }_{\lambda \nu }^\mu \right)\frac{dn^\sigma }{ds}(u^\lambda n^\rho u^\rho n^\lambda )\hfill \\ \\ +\left(\mathrm{\Gamma }_{\sigma \nu }^\tau _\tau \mathrm{\Gamma }_{\lambda \rho }^\mu +2\mathrm{\Gamma }_{\lambda \tau }^\mu _\rho \mathrm{\Gamma }_{\sigma \nu }^\tau _\nu _\sigma \mathrm{\Gamma }_{\lambda \rho }^\mu \right)(u^\lambda u^\rho n^\sigma n^\nu u^\sigma u^\nu n^\lambda n^\rho ).\hfill \end{array}$$
(16)
An equation for the 3rd-order deviation is presented in Appendix 1.
The non-manifestly covariant geodesic deviation equations are often better adapted to deriving successive approximations for geodesics close to the initial one. Starting from a given geodesic $`x^\mu (s)`$ we can solve Eq. (10) and find the first-order deviation vector $`n^\mu (s)`$. Then, inserting $`u^\mu (s)`$ and $`n^\mu (s)`$, by now completely determined, into the system (16), we can solve and find the second-order deviation vector $`b^\mu (s)`$, and subsequently for the true second-order coordinate deviation $`\delta ^2x^\mu `$, and so forth. As an example, below we describe non-circular motion, along with Kepler’s laws (in an approximate version), together with the relativistic perihelion advance, starting from a circular orbit in Schwarzschild metric.
Although for orbital motion in a Schwarzschild background we have at our disposal the exact solutions in terms of quadratures (with integrals of elliptic or Jacobi type), our approach is particularly well-suited for numerical computations, because in appropriate (Gaussian) coordinates the geodesic curves can display a very simple parametric form, and all the components of the $`4`$-velocity and other quantities reduce to constants when restricted to that geodesic.
In this case equation (10) reduces to a linear system with constant coefficients, which after diagonalization becomes a collection of harmonic oscillators, and all that remains is to find the characteristic frequencies. In the next step, we get a collection of harmonic oscillators excited by external periodic forces represented by the right-hand side of (16), which can also be solved very easily, and so forth.
In the third order, the presence of resonances giving rise to secular terms could in principle lead to instability of the orbit we started with; but this phenomenon can be dealt with by Poincaré’s method , according to which such terms can be eliminated if we admit that the frequency of the resulting solution is also slightly modified by the exterior perturbation, and can be expanded in a formal series in successive powers of the initial (small) deformation parameter.
At the end, the deviation becomes a series of powers of a small parameter containing linear combination of characteristic frequencies appearing on the right-hand side, which are entire multiples of the basic frequency, also slightly deformed. This description of planetary motion as a superposition of different harmonic motions has been first introduced by Ptolemaeus in the II century . We shall now analyse the simplest case of circular orbits in Schwarzschild geometry.

## 3 Circular orbits in Schwarzschild metric

Let us consider the geodesic deviation equation starting with a circular orbit in the field of a spherically-symmetric massive body, i.e. in the Schwarzschild metric. The circular orbits and their stability have been analyzed and studied in several papers  and books, e.g. the well-known monograph by Chandrasekhar .
The gravitational field is described by the line-element (in natural coordinates with $`c=1`$ and $`G=1`$)

$$g_{\mu \nu }dx^\mu dx^\nu =ds^2=B(r)dt^2+\frac{1}{B(r)}dr^2+r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right),$$
(17)
with

$$B(r)=1\frac{2M}{r}.$$
(18)
We recall the essential features of the solution of the geodesic equations for a test particle of mass $`m<<M`$. As the spherical symmetry guarantees conservation of angular momentum, the particle orbits are always confined to an equatorial plane, which we choose to be the plane $`\theta =\pi /2`$. The angular momentum $`J`$ is then directed along the $`z`$-axis. Denoting its magnitude per unit of mass by $`\mathrm{}=J/m`$, we have

$$\frac{d\varphi }{ds}=\frac{\mathrm{}}{r^2}.$$
(19)
In addition, as the metric is static outside the horizon $`r_+=2M`$, it allows a time-like Killing vector which guarantees the existence of a conserved world-line energy (per unit of mass $`m`$) $`\epsilon `$, such that

$$\frac{dt}{ds}=\frac{\epsilon }{1\frac{2M}{r}}.$$
(20)
Finally, the equation for the radial coordinate $`r`$ can be integrated owing to the conservation of the world-line Hamiltonian, i.e. the conservation of the absolute four-velocity:

$$\left(\frac{dr}{ds}\right)^2=\epsilon ^2\left(1\frac{2M}{r}\right)\left(1+\frac{\mathrm{}^2}{r^2}\right).$$
(21)
From this we derive a simplified expression for the radial acceleration:

$$\frac{d^2r}{ds^2}=\frac{M}{r^2}+\left(\frac{\mathrm{}^2}{r^3}\right)\left(1\frac{3M}{r}\right).$$
(22)
The equation (21) can in principle be integrated directly; indeed, the orbital function $`r(\varphi )`$ is given by an elliptic integral . However, to get directly an approximate parametric solution to the equations of motion one can also study perturbations of special simple orbits. In the following we study the problem for bound orbits by considering the first and second-order geodesic deviation equations for the special case of world lines close to circular orbits.
Observe that for circular orbits $`r=R=`$ constant, the expressions for $`dr/ds`$, Eq. (21), and $`d^2r/ds^2`$, Eq. (22), must both vanish at all times. This produces two relations between the three dynamical quantities $`(R,\epsilon ,\mathrm{})`$, showing that the circular orbits are characterized completely by specifying either the radial coordinate, or the energy, or the angular momentum of the planet. In particular, the equation for null radial velocity gives

$$\epsilon ^2=\left(1\frac{2M}{R}\right)\left(1+\frac{\mathrm{}^2}{R^2}\right).$$
(23)
Then the null radial acceleration condition (22) gives the well-known result

$$MR^2\mathrm{}^2(R3M)=0R=\frac{\mathrm{}^2}{2M}\left(1+\sqrt{1\frac{12M^2}{\mathrm{}^2}}\right),$$
(24)
leading to the requirement $`R6M`$ for stable circular orbits to exist.
With this in mind, and the explicit formulae for the Christoffel coefficients of Schwarzschild metric (given in the Appendix 3), we can establish now the four differential equations that must be satisfied by the geodesic deviation 4-vector $`n^\mu (s)`$ close to a circular orbit. We recall that on the circular orbit of radius $`R`$ (which is a geodesic in the background Schwarzschild metric) we have:

$$u^t=\frac{dt}{ds}=\frac{\epsilon }{(1\frac{2M}{R})},u^r=\frac{dr}{ds}=0,u^\varphi =\frac{d\phi }{ds}=\omega _0=\frac{\mathrm{}}{R^2},u^\theta =\frac{d\theta }{ds}=0,$$
(25)
because $`r=R=\text{const.}`$, $`\theta =\pi /2=\text{const.}`$, so that $`\mathrm{sin}\theta =1`$ and $`\mathrm{cos}\theta =0`$.

## 4 Geodesic deviation around circular orbit

It turns out that the four equations are much easier to arrive at if we use the explicit form of the first-order deviation equation (10). We get without effort the first three equations, for the components $`n^\theta `$, $`n^\varphi `$ and $`n^t`$:

$$\frac{d^2n^\theta }{ds^2}=(u^\varphi )^2n^\theta =\frac{\mathrm{}^2}{R^4}n^\theta ,$$
(26)
$$\frac{d^2n^\varphi }{ds^2}=\frac{2\mathrm{}}{R^3}\frac{dn^r}{ds},\frac{d^2n^t}{ds^2}=\frac{2M\epsilon }{R^2(1\frac{2M}{R})^2}\frac{dn^r}{ds}.$$
(27)
The deviation $`n^\theta `$ is independent of the remaining three variables $`n^t`$, $`n^r`$ and $`n^\phi `$. The harmonic oscillator equation (26) for $`n^\theta `$ displays the frequency which is equal to the frequency of the circular motion of the planet itself:

$$n^\theta (s)=n_0^\theta \mathrm{cos}(\omega _0s+\gamma )=n_0^\theta \mathrm{cos}\left(\frac{\mathrm{}}{R^2}s+\gamma \right).$$
(28)
This can be interpreted as the result of a change of the coordinate system, with a new $`z`$-axis slightly inclined with respect to the original one, so that the plane of the orbit does not coincide with the plane $`z=0`$. In this case the deviation from the plane will be described by the above solution, i.e. a trigonometric function with the period equal to the period of the planetary motion. Being a pure coordinate effect, it allows us to eliminate the variable $`n^\theta `$ by choosing $`n^\theta =0`$.
It takes a little more time to establish the equation for $`n^r`$, using Eq. (10):

$$\frac{d^2n^r}{ds^2}+2\mathrm{\Gamma }_{\lambda \rho }^ru^\lambda \frac{dn^\rho }{ds}+_\sigma \mathrm{\Gamma }_{\lambda \rho }^ru^\lambda u^\rho n^\sigma =0.$$
(29)
Taking into account that only the components $`u^t`$ and $`u^\varphi `$ of the four-velocity on the circular orbit are different from zero, and recalling that we have chosen to set $`n^\theta =0`$, too, the only non-vanishing terms in the above equation are:

$$\frac{d^2n^r}{ds^2}+2\mathrm{\Gamma }_{tt}^ru^t\frac{dn^t}{ds}+2\mathrm{\Gamma }_{\varphi \varphi }^ru^\varphi \frac{dn^\varphi }{ds}+_r\mathrm{\Gamma }_{tt}^ru^tu^tn^r+_r\mathrm{\Gamma }_{\varphi \varphi }^ru^\varphi u^\varphi n^r=0.$$
(30)
Using the identities (24) and the definitions (25), we get

$$\frac{d^2n^r}{ds^2}\frac{3\mathrm{}^2}{R^4}\left(1\frac{2M}{R}\right)n^r+\frac{2M\epsilon }{R^2}\frac{dn^t}{ds}\frac{2\mathrm{}}{R}\left(1\frac{2M}{R}\right)\frac{dn^\phi }{ds}=0.$$
(31)
The system of three remaining equations can be expressed in a matrix form:

$$\left(\begin{array}{ccc}\frac{d^2}{ds^2}& \frac{2M\epsilon }{R^2(1\frac{2M}{R})^2}\frac{d}{ds}& 0\\ \frac{2M\epsilon }{R^2}\frac{d}{ds}& \frac{d^2}{ds^2}\frac{3\mathrm{}^2}{R^4}(1\frac{2M}{R})& \frac{2\mathrm{}}{R}(1\frac{2M}{R})\frac{d}{ds}\\ 0& \frac{2\mathrm{}}{R^3}\frac{d}{ds}& \frac{d^2}{ds^2}\end{array}\right)\left(\begin{array}{c}n^t\\ n^r\\ n^\phi \end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right).$$
(32)
The characteristic equation of the above matrix is

$$\lambda ^4\left[\lambda ^2+\frac{\mathrm{}^2}{R^4}\left(1\frac{2M}{R}\right)\frac{4M\epsilon ^2}{R^4(1\frac{2M}{R})^2}\right]=0,$$
(33)
which after using the identities (23) and (24) reduces to

$$\lambda ^4\left[\lambda ^2+\frac{\mathrm{}^2}{R^4}\left(1\frac{6M}{R}\right)\right]=0,$$
(34)
so that the characteristic circular frequency is

$$\omega =\frac{\mathrm{}}{R^2}\sqrt{1\frac{6M}{R}}=\omega _0\sqrt{1\frac{6M}{R}}.$$
(35)
It is obvious that the general solution contains oscillating terms $`\mathrm{cos}(\omega s)`$ ; however, before we analyse in detail this part of solution, let us consider the terms linear in the variable $`s`$ or constants: as a matter of fact, because of the presence of first and second-order derivatives with respect to $`s`$ in the matrix operator (32), the general solution may also contain the following vector:

$$\left(\begin{array}{c}(\mathrm{\Delta }u^t)s+\mathrm{\Delta }t\\ (\mathrm{\Delta }u^r)s+\mathrm{\Delta }r\\ (\mathrm{\Delta }u^\phi )s+\mathrm{\Delta }\phi \end{array}\right).$$
(36)
When inserted into the system (32), the solution is the following: $`\mathrm{\Delta }t`$ and $`\mathrm{\Delta }\phi `$ are arbitrary; $`\mathrm{\Delta }u^r=0`$, which means that the radial velocity remains null; and

$$\frac{3\mathrm{}^2}{R^4}\left(1\frac{2M}{R}\right)\mathrm{\Delta }r=\frac{2M\epsilon }{R^2}\mathrm{\Delta }u^t\frac{2\mathrm{}}{R}\left(1\frac{2M}{R}\right)\mathrm{\Delta }u^\phi =0.$$
(37)
This condition coincides with the transformation of the initial circular geodesic of radius $`R`$ to a neighboring one, with radius $`R+\mathrm{\Delta }r`$, with the subsequent variations $`\mathrm{\Delta }u^t`$ and $`\mathrm{\Delta }u^\phi `$ added to the corresponding components of the 4-velocity in order to satisfy the condition $`g_{\mu \nu }u^\mu u^\nu =1`$ in the linear approximation. After choosing an optimal value for $`r`$, we can forget about this particular solution, as well as about the arbitrary shift in the variables $`t`$ and $`\varphi `$, and investigate the oscillating part of the solution.
We shall choose the initial phase to have (with $`n_0^r>0)`$:

$$n^r(s)=n_0^r\mathrm{cos}(\omega s).$$
(38)
What remains to be done is to compare this frequency with the fundamental circular frequency $`\omega _0=\mathrm{}/R^2`$ of the unperturbed circular orbital motion.
But this discrepancy between the two circular frequencies $`\omega `$ and $`\omega _0`$ is exactly what produces the perihelion advance, and its value coincides with the value obtained in the usual way (1) in the limit of quasi-circular orbits, i.e. when $`e^20`$: we get both the correct value and the correct sign.
Let us display the complete solution for the first-order deviation vector $`n^\mu (s)`$ which takes into account only the non-trivial degrees of freedom:

$$n^\theta =0,n^r(s)=n_0^r\mathrm{cos}(\omega s),n^\phi =n_0^\phi \mathrm{sin}(\omega s),n^t=n_0^t\mathrm{sin}(\omega s).$$
(39)
The only independent amplitude is given by $`n_0^r`$, because we have

$`n_0^t`$ $`=`$ $`{\displaystyle \frac{2M\epsilon }{R^2(1\frac{2M}{R})^2\omega }}n_0^r={\displaystyle \frac{2\sqrt{M}}{\sqrt{R}\left(1\frac{2M}{R}\right)\sqrt{1\frac{6M}{R}}}}n_0^r,`$ (40)
$`n_0^\phi `$ $`=`$ $`{\displaystyle \frac{2\mathrm{}}{R^3\omega }}n_0^r={\displaystyle \frac{2\omega _0}{R\omega }}n_0^r={\displaystyle \frac{2}{R\sqrt{1\frac{6M}{R}}}}n_0^r.`$ (41)
So the trajectory and the law of motion are given by

$`r`$ $`=`$ $`Rn_0^r\mathrm{cos}(\omega s),`$ (42)
$`\phi `$ $`=`$ $`\omega _0s+n_0^\phi \mathrm{sin}(\omega s)={\displaystyle \frac{\sqrt{M}}{R^{3/2}\sqrt{1\frac{3M}{R}}}}s+n_0^\phi \mathrm{sin}(\omega s),`$ (43)
$`t`$ $`=`$ $`{\displaystyle \frac{\epsilon }{(1\frac{2M}{R})}}s+n_0^t\mathrm{sin}(\omega s)={\displaystyle \frac{1}{\sqrt{1\frac{3M}{R}}}}s+n_0^t\mathrm{sin}(\omega s),`$ (44)
where the phase in the argument of the $`cosine`$ function was chosen so that $`s=0`$ corresponds to the perihelion, and $`s=\frac{\pi }{\omega }`$ to the aphelion. It is important to note once again that the coefficient $`n_0^r`$, which also fixes the values of the two remaining amplitudes, $`n_0^t`$ and $`n_0^\phi `$, defines the size of the actual deviation, so that the ratio $`\frac{n_0^r}{R}`$ becomes the dimensionless infinitesimal parameter controlling the approximation series with consecutive terms proportional to the consecutives powers of $`\frac{n_0^r}{R}`$.
What we see here is the approximation to an elliptic orbital movement as described by the presence of an epicycle (exactly like in the Ptolemean system , except for the fact that the Sun is placed in the center instead of the Earth). As a matter of fact, the development into power series with respect to the eccentricity $`e`$ considered as a small parameter, and truncating all the terms except the linear one, leads to the Kepler result ,

$$r(t)=\frac{a(1e^2)}{1+e\mathrm{cos}(\omega _0t)}a\left[1e\mathrm{cos}(\omega _0t)\right],$$
(45)
which looks almost as our formula (42) if we identify the eccentricity $`e`$ with $`\frac{n_0^r}{R}`$ and the greater half-axis $`a`$ with $`R`$; but there is also the additional difference, that the circular frequency of the epicycle is now slightly lower than the circular frequency of the unperturbed circular motion.
But if the circular frequency is lower, the period is slightly longer: in a linear approximation, we have

$$\omega =\sqrt{\frac{\mathrm{}^2}{R^4}\left(1\frac{6M}{R}\right)},$$
(46)
hence keeping the terms up to the third order in $`\frac{M}{R}`$,

$$TT_0\left(1+\frac{3M}{R}+\frac{27}{2}\frac{M^2}{R^2}+\frac{135}{2}\frac{M^3}{R^3}+\mathrm{}\right).$$
(47)
Then obviously one must have $`\frac{\mathrm{\Delta }\phi }{2\pi }=\frac{\mathrm{\Delta }T}{T_0}`$ from which we obtain the perihelion advance after one revolution

$$\mathrm{\Delta }\phi =\frac{6\pi M}{R}+\frac{27\pi M^2}{R^2}+\frac{135\pi M^3}{R^3}+\mathrm{}$$
(48)
It is obvious that at this order of approximation we could not keep track of the factor $`(1e^2)^1`$, containing the eccentricity (here replaced by the ratio $`\frac{n_0^r}{R}`$) only through its square. In contrast, we obtain without effort the coefficients in front of terms quadratic or cubic in $`\frac{M}{R}`$. This shows that our method can be of interest when one has to consider the low-eccentricity orbits in the vicinity of very massive and compact bodies, having a non-negligible ratio $`\frac{M}{R}`$.
In order to include this effect, at least in its approximate form as the factor $`(1+e^2)`$, we must go beyond the first-order deviation equations and investigate the solutions of the equations describing the quadratic effects (16).

## 5 The second-order geodesic deviation

After inserting the complete solution for the first-order deviation vector (39)–(41) into the system (15) and a tedious calculation, we find the following set of linear equations satisfied by the second-order deviation vector $`b^\mu (s)`$:

$$\left(\begin{array}{ccc}\frac{d^2}{ds^2}& \frac{2M\epsilon }{R^2(1\frac{2M}{R})^2}\frac{d}{ds}& 0\\ \frac{2M\epsilon }{R^2}\frac{d}{ds}& \frac{d^2}{ds^2}\frac{3\mathrm{}^2}{R^4}(1\frac{2M}{R})& \frac{2\mathrm{}}{R}(1\frac{2M}{R})\frac{d}{ds}\\ 0& \frac{2\mathrm{}}{R^3}\frac{d}{ds}& \frac{d^2}{ds^2}\end{array}\right)\left(\begin{array}{c}b^t\\ b^r\\ b^\phi \end{array}\right)=(n_0^r)^2\left(\begin{array}{c}C^t\\ C^r\\ C^\phi \end{array}\right),$$
(49)
where we have put into evidence the common factor $`(n_0^r)^2`$, which shows the explicit quadratic dependence of the second-order deviation vector $`b^\mu `$ on the first-order deviation amplitude $`n_0^r`$. The constants $`C^t,C^r`$ and $`C^\phi `$ are expressions depending on $`M`$, $`R`$, $`\omega _0`$, $`\omega `$, $`\epsilon `$, $`\mathrm{sin}(2\omega s)`$ and $`\mathrm{cos}(2\omega s)`$:

$$C^t=\frac{6M^2(2\frac{7M}{R})\epsilon \mathrm{sin}\left(2\omega s\right)}{(1\frac{3M}{R})(1\frac{2M}{R})^2R^6\omega },$$
(50)
$$C^r=\frac{3M\left[(2\frac{5M}{R}+\frac{18M^2}{R^2})(6\frac{27M}{R}+\frac{6M}{R^2}^2)\mathrm{cos}\left(2\omega s\right)\right]}{2(1\frac{3M}{R})(1\frac{6M}{R})R^4},$$
(51)
$$C^\phi =\frac{6M(1\frac{M}{R})\omega _0\mathrm{sin}\left(2\omega s\right)}{(1\frac{3M}{R})R^5\omega }.$$
(52)
The solution of the above matrix for $`b^\mu (s)`$ has the same characteristic equation of the matrix (32) for $`n^\mu (s)`$, and the general solution containing oscillating terms with angular frequency $`\omega `$ is of no interest because it is already accounted for by $`n^\mu (s)`$. But the particular solution includes the terms linear in the proper time $`s`$, constant ones, and the terms oscillating with angular frequency $`2\omega `$:

$$b^t=\frac{\left(n_0^r\right)^2M\epsilon }{R^3(1\frac{6M}{R})(1\frac{2M}{R})^2}\left[\frac{3(2\frac{5M}{R}+\frac{18M^2}{R^2})}{1\frac{6M}{R}}s+\frac{2\frac{13M}{R}}{\omega }\mathrm{sin}(2\omega s)\right],$$
(53)
$$b^r=\frac{\left(n_0^r\right)^2}{2R(1\frac{6M}{R})}\left[\frac{3(2\frac{5M}{R}+\frac{18M^2}{R^2})}{1\frac{6M}{R}}+\left(2+\frac{5M}{R}\right)\mathrm{cos}\left(2\omega s\right)\right],$$
(54)
$$b^\phi =\frac{\left(n_0^r\right)^2\omega _0}{R^2(1\frac{6M}{R})}\left[\frac{3(2\frac{5M}{R}+\frac{18M^2}{R^2})}{1\frac{6M}{R}}s+\frac{1\frac{8M}{R}}{2\omega }\mathrm{sin}(2\omega s)\right].$$
(55)
As explained in Section $`2`$, we need to calculate $`\frac{1}{2}\delta ^2x^\mu `$ to obtain the geodesic curve $`x^\mu `$ with second-order geodesic deviation:

$$\delta ^2t=\frac{\left(n_0^r\right)^2M\epsilon }{R^3}\left[\frac{3(2\frac{5M}{R}+\frac{18M}{R^2}^2)}{(1\frac{2M}{R})^2(1\frac{6M}{R})^2}s+\frac{2\frac{15M}{R}+\frac{14M^2}{R^2}}{(1\frac{6M}{R})(1\frac{2M}{R})^3\omega }\mathrm{sin}(2\omega s)\right],$$
(56)
$$\delta ^2r=\frac{\left(n_0^r\right)^2}{R(1\frac{6M}{R})}\left[\frac{5\frac{33M}{R}+\frac{90M^2}{R^2}\frac{72M^3}{R^3}}{(1\frac{2M}{R})(1\frac{6M}{R})}\left(1\frac{7M}{R}\right)\mathrm{cos}(2\omega s)\right],$$
(57)
$$\delta ^2\phi =\frac{\left(n_0^r\right)^2\omega _0}{R^2(1\frac{6M}{R})}\left[\frac{3(2\frac{5M}{R}+\frac{18M^2}{R^2})}{1\frac{6M}{R}}s+\frac{5\frac{32M}{R}}{2\omega }\mathrm{sin}(2\omega s)\right].$$
(58)
The fact that the second-order deviation vector $`b^\mu `$ turns with angular frequency $`2\omega `$ enables us to get a better approximation of the elliptic shape of the resulting orbit. The trajectory described by $`x^\mu `$ including second-order deviations is not an ellipse, but we can match the perihelion and aphelion distances to see that $`Ra`$ and $`en_0^r/R`$ when second-order deviation is used. The perihelion and aphelion distances of the Keplerian, i.e., elliptical orbit are $`a(1e)`$ and $`a(1+e)`$. For $`x^\mu `$, the perihelion is obtained when $`\omega s=2k\pi `$ and the aphelion when $`\omega s=(1+2k)\pi `$, where $`k`$. Matching the radius for perihelion and aphelion, we obtain the semimajor axis $`a`$ and the eccentricity $`e`$ of an ellipse that has the same perihelion and aphelion distances of the orbit described by $`x^\mu `$:

$$a=R+\frac{\left(n_0^r\right)^2}{12R}\left[1+\frac{3}{1\frac{2M}{R}}+\frac{7}{1\frac{6M}{R}}+\frac{15}{(1\frac{6M}{R})^2}\right]$$
(59)
$$e=\frac{n_0^r(1\frac{2M}{R})(1\frac{6M}{R})^2}{R(1\frac{2M}{R})(1\frac{6M}{R})^2+\frac{\left(n_0^r\right)^2}{R}\left[2\frac{9M}{R}+\frac{11M^2}{R^2}+\frac{6M^3}{R^3}\right]}=\frac{n_0^r}{R}+𝒪\left(\frac{(n_0^r)^3}{R^3}\right).$$
(60)
In the limit case of $`\frac{M}{R}0`$, there is no perihelion advance and $`a=R\left[1+2(\frac{n_0^r}{R})^2\right]`$ and $`e=\frac{n_0^r}{R}`$, so the second-order deviation increases the semimajor axis $`a`$ of a matching ellipse compared to the first-order deviation, when $`a=R`$ and $`e=\frac{n_0^r}{R}`$.
Another comparison with elliptic orbits concerns the shape of the orbit described by $`r(\phi )`$. From $`\phi (s)`$ it is possible to write $`s(\phi )`$ by means of successive approximations, beginning with $`\omega s=\phi \sqrt{1\frac{6M}{R}}`$. Finally, $`s`$ can be replaced in $`r(s)`$ and we obtain $`r(\phi )`$ up to the second order in $`\frac{n_0^r}{R}`$:

$`{\displaystyle \frac{r}{R}}`$ $`=`$ $`1{\displaystyle \frac{n_0^r}{R}}\mathrm{cos}\left({\displaystyle \frac{\omega }{\omega _0}}\phi \right)+\left({\displaystyle \frac{n_0^r}{R}}\right)^2[{\displaystyle \frac{3\frac{5M}{R}\frac{30M^2}{R^2}+\frac{72M^3}{R^3}}{2(1\frac{2M}{R})(1\frac{6M}{R})^2}}`$ (62)
$`+{\displaystyle \frac{(1\frac{5M}{R})}{2(1\frac{6M}{R})}}\mathrm{cos}\left({\displaystyle \frac{2\omega }{\omega _0}}\phi \right)]+\mathrm{}`$
In the limit $`\frac{M}{R}0`$, the exact equation of a ellipse is obtained up to the second order in $`e`$, where $`e=n_0^r/R`$ and $`r_0=(1+e^2)R`$:

$$r=\frac{r_0}{1+e\mathrm{cos}\phi }=\frac{(1+e^2)R}{1+e\mathrm{cos}\phi }=R\left[1e\mathrm{cos}\phi +e^2\left(\frac{3}{2}+\frac{1}{2}\mathrm{cos}2\phi \right)+\mathrm{}\right].$$
(63)
Comparing with the ellipse equation (45), we have $`r_0=a(1e^2)`$, so $`a=\frac{R(1+e^2)}{(1e^2)}R(1+2e^2)`$ which agrees with the analysis of Eqs. (59)–(60).

## 6 Third-order terms and Poincaré’s method

With the third-order approximation we are facing a new problem, arising from the presence of resonance terms on the right-hand side. It is easy to see that after reducing the expressions on the right-hand side of equation (72) in Appendix 1, which contain the terms of the form

$$\mathrm{cos}^3\omega s,\mathrm{sin}\omega s\mathrm{cos}^2\omega s$$

and the like, we shall get not only the terms containing

$$\mathrm{sin}\mathrm{\hspace{0.17em}3}\omega s,\mathrm{and}\mathrm{cos}\mathrm{\hspace{0.17em}3}\omega s$$

which do not create any particular problem, but also the resonance terms containing the functions $`\mathrm{sin}\omega s`$ and $`\mathrm{cos}\omega s`$ , whose circular frequency is the same as the eigenvalue of the matrix-operator acting on the left-hand side.

As a matter of fact, the equation for the covariant third-order deviation $`h^\mu `$ can be written in matrix form, with principal part linear in the third-order deviation $`h^\mu `$, represented by exactly the same differential operator as in the lower-order deviation equations. The right-hand side is separated into two parts, one oscillating with frequency $`\omega `$, and another with frequency $`3\omega `$:

$$\begin{array}{c}\left(\begin{array}{ccc}\frac{d^2}{ds^2}& \frac{2M\epsilon }{R^2(1\frac{2M}{R})^2}\frac{d}{ds}& 0\\ \frac{2M\epsilon }{R^2}\frac{d}{ds}& \frac{d^2}{ds^2}\frac{3\mathrm{}^2}{R^4}(1\frac{2M}{R})& \frac{2\mathrm{}}{R}(1\frac{2M}{R})\frac{d}{ds}\\ 0& \frac{2\mathrm{}}{R^3}\frac{d}{ds}& \frac{d^2}{ds^2}\end{array}\right)\left(\begin{array}{c}h^t\\ h^r\\ h^\phi \end{array}\right)=\hfill \\ \\ =(n_0^r)^3\left(\begin{array}{c}B^t\mathrm{sin}(\omega s)+C^t\mathrm{sin}(\mathrm{\hspace{0.17em}3}\omega s)+sD^t\mathrm{cos}(\omega s)\\ B^r\mathrm{cos}(\omega s)+C^r\mathrm{cos}(\mathrm{\hspace{0.17em}3}\omega s)+sD^r\mathrm{sin}(\omega s)\\ B^\phi \mathrm{sin}(\omega s)+C^\phi \mathrm{sin}(\mathrm{\hspace{0.17em}3}\omega s)+sD^\phi \mathrm{cos}(\omega s)\end{array}\right),\hfill \end{array}$$
(64)
where the coefficients $`B^k`$, $`C^k`$ and $`D^k`$, $`k=t,r,\phi `$ are complicated functions of $`\frac{M}{R}`$.

The proper frequency of the matrix operator acting on the left-hand side is equal to $`\omega `$; the terms containing the triple frequency $`3\omega `$ will give rise to the unique non-singular solution of the same frequency, but the resonance terms of the basic frequency on the right-hand side will give rise to secular terms, proportional to $`s`$, which is in contradiction with the bounded character of the deviation we have supposed from the beginning. The term proportional to $`s`$ on the right-hand side is eliminated in the differential equation for $`h^r`$ when $`\frac{dh^\phi }{ds}`$ and $`\frac{dh^t}{ds}`$ are replaced by theirs values.

Poincaré  was first to understand that in order to solve this apparent contradiction, one has to take into account possible perturbation of the basic frequency itself, which amounts to the replacement of $`\omega `$ by an infinite series in powers of the infinitesimal parameter, which in our case is the eccentricity $`e=\frac{n_0^r}{R}`$:

$$\omega \omega +e\omega _1+e^2\omega _2+e^3\omega _3+\mathrm{},$$
(65)
Then, developing both sides into a series of powers of the parameter $`e`$, we can not only recover the former differential equations for the vectors $`n^\mu ,b^\mu ,h^\mu `$, but get also some algebraic relations defining the corrections $`\omega _1,\omega _2,\omega _3,`$ etc.
The equations resulting from the requirement that all resonant terms on the right-hand side be canceled by similar terms on the left-hand side are rather complicated. We do not attempt to solve them here. However, one easily observes that the absence of resonant terms in the second-order deviation equations forces $`\omega _1`$ to vanish, while the next term $`\omega _2`$ is different from $`0.`$
Similarly, as there are no resonant terms in the equations determining the fourth-order deviation, because all four-power combinations of sine and cosine functions will produce terms oscillating with frequencies $`2\omega `$ and $`4\omega `$; as a result, the correction $`\omega _3`$ will be also equal to $`0`$. Next secular terms will appear at the fifth-order approximation, as products of the type $`cos^5\omega s`$, $`sin^3\omega scos^2\omega s`$, etc, produce resonant terms again, which will enable us to find the correction $`\omega _4`$, and so on, so that the resulting series representing the frequency $`\omega `$ contains only even powers of the small parameter $`\frac{n_0^r}{R}`$.

## 7 Gravitational radiation

The decomposition of the elliptic trajectory turning slowly around its focal point into a series of epicycles around a circular orbit can also serve for obtaining an approximate spectral decomposition of gravitational waves emitted by a celestial body moving around a very massive attracting center.
It is well known that gravitational waves are emitted when the quadrupole moment of a mass distribution is different from zero, and the amplitude of the wave is proportional to the third derivative of the quadrupole moment with respect to time (in the reference system in which the center of mass coincides with the origin of the Cartesian basis in three dimensions, see Ref. ).
Of course, it is only a linear approximation, but it takes the main features of the gravitational radiation emitted by the system well into account, provided the velocities and the gravitational fields are not relativistic and the wavelength of gravitational radiation is large compared to the dimensions of the source (quadrupole approximation).
More precisely, let us denote the tensor $`Q_{ij}`$ of a given mass distribution $`\mu (x_i),`$ where $`i,j,=1,2,3`$, see Ref. :

$$Q_{ij}=\mu x_ix_j𝑑V=\underset{\alpha }{}m_\alpha x_{\alpha i}x_{\alpha j},$$
(66)
where $`m_\alpha `$ are point masses.
Let $`\stackrel{}{OP}`$ be the vector pointing at the observer (placed at the point $`P`$), from the origin of the coordinate system coinciding with the center of mass of the two orbiting bodies whose motion is approximately described by our solution in a Fourier series form. It is also supposed that the length of this vector is much greater than the characteristic dimensions of the radiating system, i.e. $`\stackrel{}{OP}R`$.
Then the total power of gravitational radiation $`P`$ emitted by the system over all directions is given by the following expression (see Ref. ):

$$P=\frac{G}{5c^5}\left(\frac{d^3Q_{ij}}{dt^3}\frac{d^3Q_{ij}}{dt^3}\frac{1}{3}\frac{d^3Q_{ii}}{dt^3}\frac{d^3Q_{jj}}{dt^3}\right).$$
(67)
When applied to Keplerian motion of two masses $`m_1`$ and $`m_2`$, with orbit equation and angular velocity given by

$$r=\frac{a(1e^2)}{1+e\mathrm{cos}\phi },\frac{d\phi }{dt}=\frac{\sqrt{G(m_1+m_2)a(1e^2)}}{r^2},$$
(68)
the total power $`P`$ now reads

$$P=\frac{8}{15}\frac{G^4}{c^5}\frac{m_1^2m_2^2(m_1+m_2)}{a^5(1e^2)^5}(1+e\mathrm{cos}\phi )^4\left[12(1+e\mathrm{cos}\phi )^2+e^2\mathrm{sin}^2\phi \right].$$
(69)
We shall calculate the $`P`$ in Eq. (67) with our solution $`x^\mu `$ using second-order geodesic deviation, to inspect the non-negligible effects of the ratio $`\frac{M}{R}`$. We have the explicit solutions $`r(s)`$, $`\phi (s)`$ and $`t(s)`$, so to calculate $`\frac{dQ_{ij}}{dt}`$ we need only the derivatives with respect to $`s`$, i.e., $`\frac{df}{dt}=\frac{df}{ds}/\frac{dt}{ds}`$ can be applied successively to obtain $`\frac{d^3Q_{ij}}{dt^3}`$. So we finally get $`P`$ as function of $`s`$, which is not shown here because it is a very large expression that nevertheless can be easily obtained using a symbolic calculus computer program.
As we want to compare the two total powers $`P`$ during one orbital period (between perihelions), $`P`$ in the Kepler case is obtained from the numerical solution for $`\phi (t)`$ calculated from Eq. (68), and $`P`$ of the geodesic deviation case has to use $`s(t)`$ obtained from $`t(s)`$ by means of successive approximations, starting with $`s=\frac{t}{\epsilon }\sqrt{1\frac{2M}{R}}`$.
There are many possible ways to compare a Keplerian orbit with a relativistic one. Here we assume $`m_1m_2`$ and fix the values of $`a`$, $`e`$, $`m_1`$; the values of $`R`$ and $`n_0^r`$ are calculated to obtain an exact ellipse (up to the second order in $`e`$) in the limit $`\frac{M}{R}0`$, like Eq. (63), so $`R=\frac{(1e^2)}{(1+e^2)}a`$ and $`n_0^r=Re`$. Up to first order in $`e`$, we have $`R=a`$. The choice of $`M=m_1`$ allows the two total powers $`P`$ to be equal when $`e=0`$ and $`\frac{M}{R}0`$. Figures $`1`$ and $`2`$ show this comparison for small eccentricities and non-negligible $`\frac{M}{R}`$ ratios.

Because the emitted total powers $`P`$ calculated with geodesic deviations depend on the $`\frac{M}{R}`$ ratio, we see that the period is not $`T=\frac{2\pi a^{3/2}}{\sqrt{Gm_1}}`$ (third Kepler’s law), but an increased one,

$$T=\frac{2\pi R^{3/2}}{\sqrt{GM}\sqrt{1\frac{6GM}{R}}}+𝒪\left(\frac{(n_0^r)^2}{R^2}\right).$$
(70)
This effect is the direct consequence of the form of angular frequency $`\omega `$ that appears in the first and higher-order geodesic deviations.
Another expected feature of Figures $`1`$ and $`2`$ : as $`e`$ (i.e., $`\frac{n_0^r}{R}`$) is kept small, the $`P`$ using geodesic deviations converge very fast in respect of the orders of geodesic deviation.
Caution is required as the use of quadrupole approximation is not allowed for high values of $`\frac{M}{R}`$, so the exact amplitude and shape of $`P`$ using geodesic deviations can only be calculated if additional $`\frac{M}{R}`$ contributions to the gravitational radiation formula are included. This approach, but using the post-Newtonian expansion scheme, is well developed in Refs. .

## 8 Discussion

In this paper we have introduced an extension of the geodesic deviation idea in order to calculate approximate orbits of point masses in gravitational fields. This scheme is of practical applicability to the problem of the emission of gravitational radiation. Although in the present paper we restricted our investigations to the case of Schwarzschild background fields, our method can be easily extended to other background fields . An example is provided by the discussion of Reissner-Nordstrøm fields in Ref. .
Since the initial work by Einstein, the problem of orbits and radiation is addressed in the literature mostly through the post-Newtonian expansion scheme , -. In this approach the starting point for the successive approximations is found in Newtonian theory, whilst relativistic effects are introduced by corrections of higher order in $`\frac{v}{c}`$ or $`\frac{M}{R}`$. The advantage of this approach is, that one can start with an orbit of arbitrary high eccentricity. In contrast, our approach starts from a true solution of the relativistic problem, but in the case of the Schwarzschild background we choose a circular one. We then approach finite eccentricity orbits in a fully relativistic scheme, by summing up higher-order geodesic deviations, for which we have derived the explicit expressions.

The two approaches are complementary in the following sense: the post-Newtonian scheme gives better results for small values of $`\frac{M}{R}`$ and arbitrary eccentricity, whereas our scheme is best adapted for small eccentricities, but arbitrary values of $`\frac{M}{R}<\frac{1}{6}`$. In both approaches the emission of gravitational radiation is estimated using the quadrupole formula, based on a flat-space approximation.

The next challenge is to include finite-size and radiation back-reaction effects. In the post-Newtonian scheme some progress in this direction has already been made. In this aspect our result may be regarded as the first term in an expansion in $`\frac{m}{M}`$. Other applications can be found in problems of gravitational lensing and perturbations by gravitational waves.

Acknowledgments
R. K. and R. C. Jr. wish to thank Dr. Christian Klein for enlightening discussions and Cédric Leygnac for his valuable help in checking the results. We also wish to express our thanks to the Referees for their constructive criticisms and useful remarks that helped to improve the presentation.

## Appendix 1

The covariant third-order deviation equation is obtained via the same procedure that has served to derive the second-order covariant deviation. The third-order geodesic deviation itself is

$`\delta ^3x^\mu `$ $`=`$ $`(pp_0)^3{\displaystyle \frac{^3x^\mu }{p^3}}`$ (71)
$`=`$ $`ϵ^3\left[h^\mu 3\mathrm{\Gamma }_{\lambda \nu }^\mu n^\lambda b^\nu +\left(_\kappa \mathrm{\Gamma }_{\lambda \nu }^\mu 2\mathrm{\Gamma }_{\lambda \sigma }^\mu \mathrm{\Gamma }_{\kappa \nu }^\sigma \right)n^\kappa n^\lambda n^\nu \right],`$
where $`h^\mu =Db^\mu /Dp=D^2n^\mu /Dp^2`$. We derive the third-order deviation equation by taking the covariant derivative w.r.t. $`p`$ of Eq. (15) for $`b^\mu `$, with the result

$$\begin{array}{cc}\frac{D^2h^\mu }{Ds^2}\hfill & +R_{\rho \lambda \sigma }^\mu u^\lambda u^\sigma h^\rho =\mathrm{\hspace{0.17em}\; 6}R_{\lambda \rho \sigma }^\mu (u^\lambda n^\rho \frac{Db^\sigma }{Ds}+u^\lambda \frac{Dn^\sigma }{Ds}b^\rho )\hfill \\ & \\ & +_\tau _\nu R_{\lambda \rho \sigma }^\mu (u^\nu u^\rho n^\tau n^\lambda n^\sigma +u^\rho u^\sigma n^\tau n^\nu n^\lambda )\hfill \\ & \\ & +_\nu R_{\lambda \rho \sigma }^\mu [4u^\lambda n^\nu n^\rho \frac{Dn^\sigma }{Ds}+u^\nu u^\rho n^\lambda b^\sigma +u^\rho u^\sigma n^\nu b^\lambda \hfill \\ & \\ & +2(u^\rho n^\lambda n^\nu \frac{Dn^\sigma }{Ds}+u^\nu n^\lambda n^\sigma \frac{Dn^\rho }{Ds}+u^\nu u^\rho n^\lambda \frac{Db^\sigma }{Ds}+u^\rho u^\sigma n^\nu \frac{Db^\lambda }{Ds})]\hfill \\ & \\ & +4R_{\lambda \rho \sigma }^\mu n^\rho \frac{Dn^\lambda }{Ds}\frac{Dn^\sigma }{Ds}+4R_{\lambda \rho \sigma }^\mu R_{\alpha \beta \gamma }^\sigma u^\lambda u^\beta n^\gamma n^\alpha n^\rho .\hfill \end{array}$$
(72)
Related studies of higher-order differentials and their covariant generalizations from a more general perspective can be found in recent papers .

## Appendix 2

The contributions of various orders to the geodesic deviation obtained in Section $`2`$ can be deduced in an elegant, coordinate-independent and slightly more general manner . Given a one-parameter congruence of geodesics, one can define the tangent vector field $`Z`$ and the local Jacobi field $`X`$; then the Lie bracket of these fields vanishes (because the congruence spans a submanifold and therefore is integrable), so that $`[X,Z]=0`$.
The geodesic equation is $`_ZZ=0`$. Then, applying the definition of Riemann tensor to the vectors $`X,Y`$ and $`Z`$:

$$\left[_X_Z_Z_X\right]Z_{[X,Z]}Z=R(X,Z)Z,$$
(73)
and taking into account that $`[X,Z]=0`$ as well as the fact that $`_XZ=_ZX`$ and the anti-symmetry of $`R(X,Z)`$ in its two arguments, we get easily

$$_Z_XZ=_Z^2X=R(Z,X)Z,$$
(74)
which coincides with Eq. (9) for the geodesic deviation, after we identify the components of the vector fields as $`Z^\mu =u^\mu ,X^\mu =n^\mu `$.

One may continue in the same spirit and introduce two linearly independent Jacobi fields, $`X`$ and $`Y`$, both satisfying $`[X,Z]=0=[Y,Z]`$, to obtain the coordinate-independent form of Eq. (15), as follows. The two linearly independent Jacobi fields, $`X`$ and $`Y`$, satisfy $`[X,Z]=0`$ and $`[Y,Z]=0`$. By virtue of the Jacobi identity, we have $`[[X,Y],Z]=0`$, hence $`[X,Y]`$ is also a Jacobi field (i.e., it satisfies Eq. (74)). Applying the same formula to this field, we get

$$_Z^2\left([X,Y]\right)=R(Z,[X,Y])Z$$
(75)
Then, using the fact that $`_XY_YX=[X,Y]`$, we can write the left hand side of the above equation as $`_Z^2(_XY_YX])`$.

Next, we can write this equation explicitly as

$$_Z^2(_XY_YX)=R(Z,_XY_YX)Z$$

and furthermore, using the linearity property, as

$$_Z^2(_YX)R(Z,_YX)Z=_Z^2(_XY)R(Z,_XY)Z.$$
(76)
The left-hand side of the above equation coincides with the usual Jacobi equation applied to the field $`_YX`$, whereas the right-hand side can be transformed using the definition of the Riemann tensor: the term $`_Z^2(_XY)`$ gives

$$\begin{array}{ccc}_Z^2Y\hfill & =\hfill & _Z(_Z_XY)=_Z(_Z_XY_X_ZY)+_Z(_X_ZY)\hfill \\ & & \\ & =\hfill & _Z(R(Z,X)Y)+_Z(_X_ZY)\hfill \\ & & \\ & =\hfill & \left[(_ZR)(Z,X)\right]Y+R(Z,_ZX)Y+R(Z,X)_ZY+_Z(_X_ZY)\hfill \end{array}$$
(77)
(here we used the fact that $`_ZZ=0`$). Manipulating further in the same manner the commutators of covariant derivations, we arrive at the result

$$_Z^2(_YX)R(Z,_YX)Z=$$

$$_XR(Z,Y)Z+_ZR(Z,X)Y+2R(Z,Y)_ZX+2R(Z,X)_ZY$$
(78)
which is equivalent to Eq. (15) upon identification $`Y=X`$ and $`b^\mu =(_XX)^\mu `$.

Although these coordinate-independent derivations are more elegant, their results are not so useful for pratical computations, i.e., the non-manisfestly covariant form of the results is better adapted for the calculus of successive deviations in a given local coordinate system.

## Appendix 3

Connections and curvatures for Schwarzschild geometry
In this appendix we collect the expressions for the components of the connections and Riemann curvature used in the main body of the paper. A. Connections. From the line-element (17) one derives the following expressions for the connection coefficients:

$$\mathrm{\Gamma }_{\rho \lambda }^\mu =\frac{1}{2}g^{\mu \sigma }(_\rho g_{\sigma \lambda }+_\lambda g_{\rho \sigma }_\sigma g_{\rho \lambda });$$
(79)
$$\begin{array}{ccc}\mathrm{\Gamma }_{rt}^t\hfill & =\hfill & \mathrm{\Gamma }_{rr}^r=\frac{M}{r^2\left(1\frac{2M}{r}\right)},\mathrm{\Gamma }_{tt}^r=\frac{M}{r^2}\left(1\frac{2M}{r}\right),\hfill \\ & & \\ \mathrm{\Gamma }_{\phi \phi }^r\hfill & =\hfill & r\mathrm{sin}^2\theta \left(1\frac{2M}{r}\right),\mathrm{\Gamma }_{\theta \theta }^r=r\left(1\frac{2M}{r}\right),\hfill \\ & & \\ \mathrm{\Gamma }_{r\theta }^\theta \hfill & =\hfill & \mathrm{\Gamma }_{r\phi }^\phi =\frac{1}{r},\mathrm{\Gamma }_{\theta \phi }^\phi =\frac{\mathrm{cos}\theta }{\mathrm{sin}\theta },\mathrm{\Gamma }_{\phi \phi }^\theta =\mathrm{sin}\theta \mathrm{cos}\theta .\hfill \end{array}$$
(80)
B. Curvature components. The corresponding curvature two-form components $`R_{\mu \nu }=\frac{1}{2}R_{\kappa \lambda \mu \nu }dx^\kappa dx^\lambda `$ are:

$$\begin{array}{ccc}R_{tr}\hfill & =\hfill & \frac{2M}{r^3}dtdr,R_{t\theta }=\frac{M}{r}\left(1\frac{2M}{r}\right)dtd\theta ,\hfill \\ & & \\ R_{t\phi }\hfill & =\hfill & \frac{M}{r}\left(1\frac{2M}{r}\right)\mathrm{sin}^2\theta dtd\phi ,R_{r\theta }=\frac{M}{r\left(1\frac{2M}{r}\right)}drd\theta ,\hfill \\ & & \\ R_{r\phi }\hfill & =\hfill & \frac{M}{r(1\frac{2M}{r})}\mathrm{sin}^2\theta drd\phi ,R_{\theta \phi }=(2Mr)\mathrm{sin}^2\theta d\theta d\phi .\hfill \end{array}$$
(81)


format: text

text:
# Cepheid and SNIa Distance Scales.

## 1. Status of Extragalactic Distance Scale

One major motivation for studying Hubble’s Constant is the complicated nature of the current standard model in cosmology, $`\mathrm{\Lambda }`$-CDM. In this model, to order of magnitude, $`\mathrm{\Omega }_{baryon}\mathrm{\Omega }_{CDM}\mathrm{\Omega }_\mathrm{\Lambda }`$ and this seems unnatural. The coincidence between the CDM and Baryon densities worried a few when CDM was first introduced (Peebles 1984, Shanks 1985). The coincidence between $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and the others worried many more (eg Dolgov, 1983, Peebles and Ratra, 1988 and Wetterich, 1988). These fine-tuning problems of the standard model are compounded by the fact that the inflation model on which the standard model sits, was partly based on a fine-tuning argument, the flatness-problem; to begin by eliminating one fine tuning problem only to end up with several gives the appearance, at least, of circular reasoning!

Shanks (1985, 1991, 1999) noted that a simpler model immediately became available if H0 actually lay below 50 kms-1 Mpc-1. An inflationary model with $`\mathrm{\Omega }_{baryon}`$=1 is then better placed to escape the baryon nucleosynthesis constraint. Simultaneously, the low value of H0 means that the X-ray gas in the Coma cluster increases towards the Coma virial mass and the lifetime of an Einstein-de Sitter Universe extends to become compatible with the ages of the oldest stars. Given the historical uncertainty there has been in the value of H0, this provides clear motivation for investigating the distance scale route to a better determination of Hubble’s Constant.

## 2. A New Era for Determining H0

Some 25 galaxies have had Cepheids detected by HST. Seventeen of these were observed by the HST Distance Scale Key Project (Freedman et al, 1994). Seven were observed in galaxies with SNIa by Sandage and collaborators (eg Sandage et al, 1996) and M96 in the Leo I Group was observed by Tanvir et al (1995). In Fig. 1 we use these data to update the comparison of I-band TF distances of Pierce & Tully (1992) with HST Cepheid distances. As can be seen, the result implies that TF distance moduli at Virgo underestimated by $``$25$`\pm `$5%. This reduces Tully-Fisher estimates of H0 from $``$85 to $``$65kms-1Mpc-1 (Giovanelli et al, 1997, Shanks 1997, Shanks, 1999, Sakai et al, 1999). The correlation of Cepheid residuals with line-width suggests TF distances may be Malmquist biased - possibly implying a bigger TF scale error at larger distances. This clear problem for TF distances, which previously has been the ‘gold standard’ of secondary distance indicators, warns that errors in the extragalactic distance scale may still be seriously underestimated!

## 3. NGC7790 Cepheid metallicity dependence?

New JKT 1.0m + CTIO 0.9m + UKIRT UBVK photometry of Cepheid Open Clusters by Hoyle et al. (2001) has uncovered an anomaly in the NGC7790 UBV 2-colour diagram, in that the F stars in the cluster show a strong UV excess with respect to zero-age main sequence stars (see Fig. 2). The result is confirmed by independent photometric data (Fry, 1997, Fry and Carney, 1997) as shown in Fig. 7b of Hoyle et al (2001). If the UV excess is caused by metallicity then NGC7790 would have \[Fe/H\]$``$ -1.5 ! To keep the Galactic Cepheid P-L relation as tight as previously observed implies that Cepheids may have a stronger metallicity dependence, $`\mathrm{\Delta }M0.66\mathrm{\Delta }[Fe/H]`$, than previously expected, in the sense that low metallicity Cepheids are intrinsically fainter. Currently we are obtaining metallicities for the F stars in NGC7790 in order to confirm this result.

## 4. HST Cepheid metallicity dependence

Meanwhile, Allen & Shanks (2001) have found an $`3\sigma `$ correlation between dispersion around the Cepheid P-L relation and galaxy metallicity for HST Cepheid galaxies (see Fig. 3). This again suggests a strong Cepheid P-L metallicity dependence and tends to support the results described above. The more extended star-formation history of high metallicity galaxies may leave a wide range of metallicities than in low metallicity galaxies like the LMC, resulting in a higher P-L dispersion if Cepheid luminosities at given period depend strongly on metallicity.

Allen & Shanks (2001) also obtain Cepheid distances via truncated maximum likelihood P-L fits to account for magnitude incompleteness caused by the non-negligible scatter in the HST P-L relations. They found that Cepheid galaxy distances at the limit of HST reach are too low. The higher than expected P-L dispersion for distant, metal-rich galaxies accentuates this effect. The conclusion is that current HST Cepheid distance moduli may be underestimated by more than 0.5 mag at the redshift of Virgo and Fornax due to both metallicity and statistical incompleteness bias. The TF distances discussed above are then underestimates by approximately 1 magnitude.

Eight HST Cepheid galaxies also have Type Ia distances. Correcting the Cepheid scale for metallicity and incompleteness bias as above and then using these distances to derive peak luminosities using the SNIa data from Gibson et al (2000), implies a strong correlation between Type Ia peak luminosity and metallicity. Such a scatter in SNIa luminosities could easily be disguised by magnitude selection effects at moderate redshifts. At higher redshift the correlation is in the right direction to explain away the need for a cosmological constant in the Supernova Hubble Diagram results, since galaxies at high redshift might be expected to have lower metallicity. Thus the conclusion is that if Cepheids have strong metallicity dependence then so have SNIa and therefore SNIa estimates of q0 and H0 may require significant corrections for metallicity.

## 5. Conclusions - Implications for H0 and SNIa

Our conclusions are as follows:-

* Key Project HST Cepheid distances imply Tully-Fisher distances at Virgo/Fornax are underestimated by $`25\pm 5`$%, reducing H0 from $``$85 to $``$65kms-1Mpc-1.
* TF distances may be Malmquist biased, suggesting there may be a bigger TF scale error at larger distances.
* If the UV excess of F stars in open cluster NGC7790 is caused by low metallicity then Cepheids have a strong metallicity dependence, $`\mathrm{\Delta }M0.66\mathrm{\Delta }[Fe/H]`$.
* Current HST Cepheid distances may be significantly underestimated at Virgo/Fornax redshifts due to metallicity and magnitude incompleteness bias, implying that values of H$`{}_{0}{}^{}<`$50kms-1Mpc-1 may still not be ruled out.
* If Cepheids have a strong metallicity dependence then so have SNIa . Thus significant metallicity corrections may need to be applied to the Type Ia Hubble Diagram before reliable estimates of q0 or H0 can be made.

### Acknowledgments.

We thank the HST Distance Scale Key Project for making the Cepheid data freely available.

## References

Allen, P.D. & Shanks, T. 2001, MNRAS, submitted.

Dolgov, A.D., 1983, In ‘The Very Early Universe’, Eds. Gibbons, G.W., Hawking, S., & Siklos, S.T.C. CUP, pp. 449-458

Freedman, W. L. et al., 1994, ApJ, 427, 628.

Fry, A.M., 1997, PhD thesis, Univ. of North Carolina.

Fry, A.M. & Carney, B.W., 1997, AJ, 113, 1073.

Gibson, B. K. et al, 1999, ApJ, 512, 48.

Giovanelli, R., Haynes, M.P., da Costa, L.N., Freudling, W., Salzer, J.J., Wegner, G., 1997, ApJ, 477, L1.

Hoyle, F., Shanks, T. & Tanvir, N.R., 2001 MNRAS submitted, astro-ph/0002521.

Peebles, P.J.E., 1984, ApJ, 284, 439.

Peebles, P.J.E. & Ratra, B. 1988, ApJ, 325, L17.

Pierce, M.J. and Tully, R.B., 1992 ApJ, 387, 47.

Sakai, S. et al, 1999, ApJ, 523, 540.

Sandage, A.R., Saha, A., Tammann, G. A., Labhardt, L., Panagia, N., Macchetto, F. D., ApJ, 1996, 460, L15.

Shanks, T., 1985, Vistas in Astronomy, 28, 595.

Shanks, T. et al., 1991, In ‘Observational Tests of Cosmological Inflation’, Eds. Shanks, T., Banday, A.J., Ellis, R.S., Frenk, C.S. & Wolfendale, A.W. pp. 205-210 Dordrecht:Kluwer.

Shanks, T., 1997, MNRAS, 290, L77.

Shanks, T., 1999, In Harmonising Cosmic Distance Scales, eds Egret D. & Heck, A., PASP, pp. 230-234.

Tanvir, N. R., Shanks, T., Ferguson, H.C. & Robinson, D.R.T. 1995, Nature, 377, 27.

Wetterich, C., 1988, Nucl. Phys. B302, 645.


text:
# Spin-related transport in ultra small Si quantum dots

## 1 Introduction

Recent advances in nanotechnology has made it possible to fabricate quantum dots so small that Coulomb blockade can be observed up to room temperature. In such small dots, single–particle energies due to the size quantization are comparable to the electrostatic charging energy and interaction effects can modify the entire energy spectrum. Important information about interactions in such quantum systems can be obtained by studying the spin states. However, it proved to be a formidable task to measure the spin of a few-electron system experimentally. So far, the most successful mapping of spins in a few-electron dot has been achieved indirectly, by comparing an experimentally obtained addition spectrum to the theoretically calculated energy spectrum for a particular geometryashoori93 ; tarucha96 . The problem of measuring spin directly in the most versatile and well studied vertical and lateral quantum dots is that the magnetic field $`B`$ dependence of their energy levels is dominated by orbital effects due to the weak confinement of electrons by electrostatic gates.

In our small Si samples we are able to measure the spin of a few-electron quantum dot directly. The dots uniquely combine the flexibility to change the number of electrons $`N`$ starting from 0 with the strong electron confinement provided by the sharp Si/SiO2 interface. At $`B<13`$ T, the $`B`$-dependence of energy levels is dominated by the Zeeman shift and we are able to measure the difference between the spin of the successive ground states $`\mathrm{\Delta }S=S(N)S(N1)`$ directly as electrons are added into the dot one by one. Also, we can follow the change of the spin configuration of the dot with a fixed number of electrons as a function of $`B`$. For one of the samples we can explain observed transition for the groundstates with up to 5 electrons within a simple model of a few energy levels which cross each other in magnetic filed as a result of Zeeman splitting. Some many-body effects, such as spontaneous dot polarization and spin blockadeweinman95  due to the spin change $`\mathrm{\Delta }S>1/2`$ has been observed.

## 2 Samples

The measurements were performed on small Si quantum dot fabricated from a silicon-on-insulator wafer. The dot resides inside a narrow bridge patterned from the top Si layer. A 50 nm thick layer of thermal oxide is grown around the bridge followed by a poly-Si gate. The fabrication steps have been described previously in detailsleobandung95 . Spacing between excited levels $`\delta 0.54`$ meV, measured using non-zero bias spectroscopy, is comparable to the charging energy $`U_c=e/C10`$ meV and is consistent with the lithographical size of the dot $`l\sqrt{\mathrm{}/m^{}\delta }100190`$ Å.

## 3 Results and discussion

Peak position as a function of $`V_g`$ is determined by the degeneracy condition that the electrochemical potentials for the ground states with $`N1`$ and $`N`$ electrons in the dot are equal. It has long been realized that for non-interacting electrons the field dependence of the peak positions $`V_g^p(B)`$ can be directly mapped onto the single-particle energy spectrum of the dot $`E(N,B)`$, provided that the Fermi energy $`E_F`$ in the contacts is field independent. For a dot with a weak confining potential we expect energy levels to shift by $`\mathrm{}\omega _c`$ or by the level spacing, whichever energy is smaller, and the shift should strongly depend on the direction of the magnetic field due to the dot anisotropy.

In Fig. 1a the conductance is plotted as a function of the gate voltage for sample E5-5b. The dot contains a few electrons (we estimate $`15`$ electrons at $`V_g=1`$ V) and the number of electrons can be changed by varying the gate voltage, each peak corresponding to an electron added to or removed from the dot. In Fig. 1b, the shift of each peak at $`B=10`$ T, relative to its zero-field position, is plotted. The average peak shift is much less than $`\mathrm{}\omega _c=6`$ meV ($`\mathrm{\Delta }V_g=52`$ mV for this sample) and is determined by the level spacing. The most striking result is that the shift is almost independent on the $`B`$ direction, which suggests that spin effects dominate.

It is possible to separate orbital and spin effects if contacts are spin polarized, as shown in Fig. 2. In the plot an excitation spectrum of another dot, E5-3b, is shown for the entrance of the first electron. At low gate voltage $`V_g<0.1`$ V the contact is spin polarized for $`B>1`$ T and the Zeeman shift of the first energy level is fully compensated by the Zeeman shift of the Fermi energy in the contacts. Thus, the resulting shift is solely due to orbital effects (assuming that $`g`$-factors in the dot and in the contacts are the same). The net shift of the levels in a parallel field (the smallest cross section) is negligible, while there is a small but measurable parabolic dependence on $`B_{}`$ (the largest cross section). The net shift at 10 T is $`<0.4`$ meV and Zeeman energy should dominate the $`B`$ dependence ($`\frac{1}{2}g^{}\mu _BB=0.6`$ meV at $`B=10`$ T). It is interesting to notice, that although orbital effects are not effecting the position of the energy levels appreciably, they can considerably change transmission coefficients. For example, at high perpendicular field the peaks, marked with the dot and triangle in Fig. 2, become suppressed, while the third peak is enhanced, compared to the $`B=0`$ value.

For a larger number of electrons $`N>3`$ the contacts are spin-degenerate up to the highest experimental field of 13 T and the shift of energy levels mainly reflects the Zeeman shift of the tunneling electron. Peaks shift with a slope $`\frac{1}{2}g^{}\mu _B`$ as a function of $`B`$ and the slope changes sign every time the two levels with different spin cross each other. Evolution of the peaks for up to 5 electrons in the dot can be understood within a model of non-interacting electrons and singlet-triplet and triplet-polarized transitions are identified as a function of $`B`$. Peak 6 has a large slope $`\frac{3}{2}g^{}\mu _B`$ at low $`B`$ and is strongly suppressed at low $`T`$. We conclude that the $`N=6`$ ground state is spontaneously polarized at $`B=0`$ and the peak suppression is due to the long predicted spin blockadeweinman95 . For larger numbers of electrons, $`N>20`$, we observed transitions with $`\mathrm{\Delta }S=3/2`$ with no apparent suppression of the peaks at low $`T`$. The spin scattering mechanism involved is not currently understood. These data will be published elsewhererokhinson00b .

The work was supported by ARO, ONR and DARPA.