You may want to use perturbation theory. This only gives you an approximate answer, but allows for analytic treatment. Your force is considered a small perturbation to the Keplerian elliptic orbit and the resulting equations of motion are expanded in powers of $K$. For linear perturbation theory, only terms linear in $K$ are retained. This simply leads to integrating the perturbation along the unperturbed original orbit. Writing your force as a vector, the perturbing acceleration is
$$
boldsymbol{a} = K frac{GM}{r^2c^2}v_rboldsymbol{v}_t
$$
with $v_r=boldsymbol{v}{cdot}hat{boldsymbol{r}}$ the radial velocity
($boldsymbol{v}equivdot{boldsymbol{r}}$) and
$boldsymbol{v}_t=(boldsymbol{v}-hat{boldsymbol{r}}(boldsymbol{v}{cdot}hat{boldsymbol{r}}))$ the rotational component of velocity (the full velocity minus the radial velocity). Here, the dot above denotes a time derivative and a hat the unit vector.
Now, it depends what you mean with 'effect'. Let's work out the changes of the orbital semimajor axis $a$, eccentricity $e$, and direction of periapse.
To summarise the results below: semi-major axis and eccentricity are unchanged, but the direction of periapse rotates in the plane of the orbit at rate
$$
omega=Omega frac{v_c^2}{c^2}
frac{K}{1-e^2},
$$
where $Omega$ is the orbital frequency and $v_c=Omega a$ with $a$ the semi-major axis. Note that (for $K=3$) this agrees with the general relativity (GR) precession rate at order $v_c^2/c^2$ (given by Einstein 1915 but not mentioned in the original question).
change of semimajor axis
From the relation $a=-GM/2E$ (with $E=frac{1}{2}boldsymbol{v}^2-GMr^{-1}$ the orbital energy) we have for the change of $a$ due to an external (non-Keplerian) acceleration
$$
dot{a}=frac{2a^2}{GM}boldsymbol{v}{cdot}boldsymbol{a}.
$$
Inserting $boldsymbol{a}$ (note that $boldsymbol{v}{cdot}boldsymbol{v}_t=h^2/r^2$ with angular momentum vector $boldsymbol{h}equivboldsymbol{r}wedgeboldsymbol{v}$), we get
$$
dot{a}=frac{2a^2Kh^2}{c^2}frac{v_r}{r^4}.
$$
Since the orbit average $langle v_r f(r)rangle=0$ for any function $f$ (see below), $langledot{a}rangle=0$.
change of eccentricity
From $boldsymbol{h}^2=(1-e^2)GMa$, we find
$$
edot{e}=-frac{boldsymbol{h}{cdot}dot{boldsymbol{h}}}{GMa}+frac{h^2dot{a}}{2GMa^2}.
$$
We already know that $langledot{a}rangle=0$, so only need to consider the first term. Thus,
$$
edot{e}=-frac{(boldsymbol{r}wedgeboldsymbol{v}){cdot}(boldsymbol{r}wedgeboldsymbol{a})}{GMa}
=-frac{r^2;boldsymbol{v}{cdot}boldsymbol{a}}{GMa}
=-frac{Kh^2}{ac^2}frac{v_r}{r^2},
$$
where I have used the identity
$(boldsymbol{a}wedgeboldsymbol{b}){cdot}(boldsymbol{c}wedgeboldsymbol{d})
=boldsymbol{a}{cdot}boldsymbol{c};boldsymbol{b}{cdot}boldsymbol{d}-
boldsymbol{a}{cdot}boldsymbol{d};boldsymbol{b}{cdot}boldsymbol{c}$ and the fact $boldsymbol{r}{cdot}boldsymbol{a}_p=0$.
Again $langle v_r/r^2rangle=0$ and hence $langledot{e}rangle=0$.
change of the direction of periapse
The eccentricity vector
$
boldsymbol{e}equivboldsymbol{v}wedgeboldsymbol{h}/GM - hat{boldsymbol{r}}
$
points (from the centre of gravity) in the direction of periapse, has magnitude $e$, and is conserved under the Keplerian motion (validate all that as an exercise!). From this definition we find its instantaneous change due to external acceleration
$$
dot{boldsymbol{e}}=
frac{boldsymbol{a}wedge(boldsymbol{r}wedgeboldsymbol{v})
+boldsymbol{v}wedge(boldsymbol{r}wedgeboldsymbol{a})}{GM}
=frac{2(boldsymbol{v}{cdot}boldsymbol{a})boldsymbol{r}
-(boldsymbol{r}{cdot}boldsymbol{v})boldsymbol{a}}{GM}
=frac{2K}{c^2}frac{h^2v_rboldsymbol{r}}{r^4}
-frac{K}{c^2}frac{v_r^2boldsymbol{v}_t}{r}
$$
where I have used the identity
$boldsymbol{a}wedge(boldsymbol{b}wedgeboldsymbol{c})=(boldsymbol{a}{cdot}boldsymbol{c})boldsymbol{b}-(boldsymbol{a}{cdot}boldsymbol{b})boldsymbol{c}$
and the fact $boldsymbol{r}{cdot}boldsymbol{a}=0$. The orbit averages of these expression are considered in the appendix below. If we finally put everything together, we get
$
dot{boldsymbol{e}}=boldsymbol{omega}wedgeboldsymbol{e}
$
with [corrected again]
$$
boldsymbol{omega}=Omega K frac{v_c^2}{c^2}
(1-e^2)^{-1}, hat{boldsymbol{h}}.
$$
This is a rotation of periapse in the plane of the orbit with angular frequency $omega=|boldsymbol{omega}|$. In particular $langle edot{e}rangle=langleboldsymbol{e}{cdot}dot{boldsymbol{e}}rangle=0$ in agreement with our previous finding.
Don't forget that due to our usage of first-order perturbation theory these results are only strictly true in the limit $K(v_c/c)^2to0$. At second-order perturbation theory, however, both $a$ and/or $e$ may change. In your numerical experiments, you should find that the orbit-averaged changes of $a$ and $e$ are either zero or scale stronger than linear with perturbation amplitude $K$.
disclaimer No guarantee that the algebra is correct. Check it!
Appendix: orbit averages
Orbit averages of $v_rf(r)$ with an abitrary (but integrable) function $f(r)$ can be
directly calculated for any type of periodic orbit. Let $F(r)$ be the antiderivative of $f(r)$, i.e. $F'!=f$, then the orbit average is:
$$
langle v_r f(r)rangle = frac{1}{T}int_0^T v_r(t),f!left(r(t)right) mathrm{d}t
= frac{1}{T} left[Fleft(r(t)right)right]_0^T = 0
$$
with $T$ the orbital period.
For the orbit averages required in $langledot{boldsymbol{e}}rangle$, we must dig a bit deeper. For a Keplerian elliptic orbit
$$
boldsymbol{r}=aleft((coseta-e)hat{boldsymbol{e}}+sqrt{1-e^2}sineta,hat{boldsymbol{k}}right)qquadtext{and}qquad
r=a(1-ecoseta)
$$
with eccentricity vector $boldsymbol{e}$ and $hat{boldsymbol{k}}equivhat{boldsymbol{h}}wedgehat{boldsymbol{e}}$ a vector perpendicular to $boldsymbol{e}$ and $boldsymbol{h}$. Here, $eta$ is the eccentric anomaly, which is related to the mean anomaly $ell$ via
$
ell=eta-esineta,
$
such that $mathrm{d}ell=(1-ecoseta)mathrm{d}eta$ and an orbit average becomes
$$
langlecdotrangle = (2pi)^{-1}int_0^{2pi}cdot;mathrm{d}ell = (2pi)^{-1}int_0^{2pi}cdot;(1-ecoseta)mathrm{d}eta.
$$
Taking the time derivative (note that $dot{ell}=Omega=sqrt{GM/a^3}$ the orbital frequency) of $boldsymbol{r}$, we find for the instantaneous (unperturbed) orbital velocity
$$
boldsymbol{v}=v_cfrac{sqrt{1-e^2}coseta,hat{boldsymbol{k}}-sineta,hat{boldsymbol{e}}}{1-ecoseta}
$$
where I have introduced $v_cequivOmega a=sqrt{GM/a}$, the speed of the circular orbit with semimajor axis $a$. From this, we find the radial velocity $v_r=hat{boldsymbol{r}}{cdot}boldsymbol{v}=v_c esineta(1-ecoseta)^{-1}$
and the rotational velocity
$$
boldsymbol{v}_t = v_cfrac{sqrt{1-e^2}(coseta-e),hat{boldsymbol{k}}-(1-e^2)sineta,hat{boldsymbol{e}}}{(1-ecoseta)^2}.
$$
With these, we have [corrected again]
$$
leftlangle frac{h^2v_rboldsymbol{r}}{r^4}rightrangle =
Omega v_c^2,hat{boldsymbol{k}},
frac{e(1-e^2)^{3/2}}{2pi}int_0^{2pi}frac{sin^2!eta}{(1-ecoseta)^4}mathrm{d}eta
=frac{Omega v_c^2e}{2(1-e^2)}hat{boldsymbol{k}}
\
leftlangle frac{v_r^2boldsymbol{v}_t}{r}rightrangle = Omega v_c^2,
hat{boldsymbol{k}},
frac{e^2(1-e^2)^{1/2}}{2pi}int_0^{2pi}frac{sin^2!eta(coseta-e)}{(1-ecoseta)^4}mathrm{d}eta=0,
$$
in particular, the components in direction $hat{boldsymbol{e}}$ average to zero. Thus [corrected again]
$$leftlangle 2frac{h^2v_rboldsymbol{r}}{r^4}-frac{v_r^2boldsymbol{v}_t}{r}rightrangle
=frac{Omega v_c^2e,hat{boldsymbol{k}}}{(1-e^2)}
$$
No comments:
Post a Comment