How does adding a inverse cube force field to a inverse square force field change the orbit?

Question

Let us say that we modify the law of gravitation to the following force field

$$ \vec{F}(r) = \frac{-GMm}{r^2}\hat{r} + \alpha\frac{|L|^2}{m^2 r^3}\hat{r} $$

where,

$L = m_2(\vec{r} \times \dot{\vec{r}} )$ which is the angular momentum of $m$ in the reference frame of $M$.

$M$ is mass of parent body which is stationery in our reference frame.

$m$ is mass of the moving body.

$\vec{r}$ is position vector of $m$ with $M$ as the reference (origin).

$\dot{\vec{r}}$ is the velocity of $m$ with M being the reference or origin again.

---

First Question

Find the equation of orbit of a planet around the Sun and draw an approximate orbit for $\alpha \ll 1.$

---

Second Question

For the Sun-Mercury system find the value of $\alpha$. Here $M$ is the Sun.

---

I could not solve this question. However I feel that Newton's Theorem of revolving Orbits is somehow needed for the solution. I would like to see a proof of the theorem as well as a complete solution for the problem.

I did not post this question for Physics SE mainly because (i) this question mainly needs knowledge of vector calculus (or so I feel) and (ii) My experience is that Math SE is a much better community so I prefer posting a borderline Maths-Physics question on Mathematics SE rather than Physics SE.

---

Edit 1

I know how to derive the two body problem by using the Rung Lenz vector method. So any answers may skip over details such as how to prove $\vec{L}$ is constant etc.

---

Edit 2

An article by Baez on the inverse cube force law.

---

Major Edit/Typo

I had made a mistake while writing the title. I had mistakenly written inverse cube potential instead of inverse cube force field.

Answer 1

Let's start by writing the kinetic energy for the planet (strictly speaking you also need to consider the motion of the sun, but for mercury $M_{\rm sun} \gg M_{\rm planet}$ so I will ignore the motion of the center of mass)

$$ T(r,\theta,\phi) = \frac{1}{2}m(\dot{r} + r^2\dot{\theta}^2 + r^2\sin^2\theta\dot{\phi}^2) \tag{1} $$

And the potential energy is just

$$ V(r) = -G\frac{Mm}{r^2}+ \alpha\frac{l^2}{2m^2 r^2} \tag{2} $$

Note that ${\bf F} = -\nabla V$, that's a good check. The Lagrangian can then be written as (down below you'll see this is completely equivalent to writing Newton's equations)

$$ L = T - V = \frac{1}{2}m(\dot{r}^2 + r^2\dot{\theta}^2 + r^2\sin^2\theta\dot{\phi}^2) + G\frac{Mm}{r^2} - \alpha\frac{l^2}{2m^2 r^2} \tag{3} $$

Now we just need to derive the equations of motion. In general if $q$ is a generalized coordinate ($q = \{r,\theta,\phi\}$) then

$$ \frac{{\rm d}}{{\rm d}t}\left(\frac{\partial L}{\partial \dot{q}}\right) - \frac{\partial L}{\partial q} = 0 \tag{4} $$

which is known as the Euler-Lagrange equation $^\color{red}{1}$.

$q = \phi$

$$ \frac{{\rm d}}{{\rm d}t}\left(\frac{\partial L}{ \partial \dot{\phi}}\right) - \frac{\partial L}{\partial \phi} = \frac{{\rm d}}{{\rm dt}}(mr^2\sin^2\theta \dot{\phi}) - (0) = 0 $$

From this you conclude that the term in parenthesis is a constant, which is nothing more than the angular momentum $l$. I will not go into the details (this is getting pretty long) but this ensures that the motion will be on a plane: $\theta = \pi/2$, we have then

$$ \frac{{\rm d}}{{\rm d}t}(mr^2\dot{\phi}) = 0 ~~~\Rightarrow~~~ l = mr^2\dot{\phi} \tag{5} $$

$q = r$

For this case we get (remember to consider the fact that $l$ is a constant)

$$ \frac{{\rm d}}{{\rm d}t}(m\dot{r})-\left(-G\frac{M m}{r^2} + \alpha\frac{l^2}{m^2r^3}+ mr \dot{\phi}^2\right) = 0 $$

Which reduces to

$$ m\frac{{\rm d}}{{\rm d}t}\dot{r} = - G\frac{mM}{r^2} + \frac{l^2}{mr^3}\left(1 - \frac{\alpha}{m} \right) \tag{6} $$

You can gain some insight by realizing this can be written as

$$ m\frac{{\rm d}}{{\rm d}t} \dot{r} = -\frac{{\rm d}}{{\rm d}r}V_{\rm eff}(r) $$

for some potential $V_{\rm eff}$ (please derive it) and finding the stable points. The variable $r$ will then more around such point

I will leave the remaining details for you to complete. Good luck!

---

$^\color{red}{1}$ You can convince yourself that this reduces to Newton's equations with a simple example, in 1D $T = m\dot{x}^2/2$ so that $L_{1D} = m\dot{x}^2/2-V$ and Eq. (4) becomes

$$ \frac{{\rm d}}{{\rm d}t}(m \dot{x}) - (-V'(x)) = 0 ~~~\Rightarrow~~~ m\ddot{x} = -V'(x) = F(x) $$

Answer 2

This is a problem with two bodies attracting each other by a conservative central force. You can treat this problem by considering the movement of just one body of reduced mass $\mu = \frac{m_1 m_2}{m_1 + m_2}$ moving in a plane - in this case, it is not difficult to show that the angular momentum of the system is preserved, and therefore the movement happens in a plane - with position vector $\vec{r} = x\vec{i} + y\vec{j}$. You can think of this movement as if the referential was one of the particles and the other had mass $\mu$, and the position vector is their relative position. To find the parametrization of the trajectory, in polar coordinates of the plane where the movement occurs, solve for $u = 1/r$ the following equation: $$\frac{d^2 u}{d \theta} + u = - \frac{\mu}{l^2}\frac{1}{u^2}F(1/u)$$ Here, $l$ denotes the magnitude of the angular momentum (which I'm assuming to be constant in this case).

For sources on this, you have Goldstein's excellent "Classical Mechanics" and Marion's and Thornton's "Classical Dynamics of Particles and Systems" (chapters 3 and 8, respectively). Now, denoting your force by $F(r) = -\frac{k_1}{r^2} - \frac{k_2}{r^3}$ (the vector force is $\vec{F}(r) = F(r) \frac{\vec{r}}{r}$), the differential equation becomes: $$\frac{d^2 u}{d \theta} + \alpha u = \beta$$ Where $\alpha$ and $\beta$ are constants in function of $\mu$, $l$, $k_1$ and $k_2$.

If $\alpha > 0$, the solution has the form: $\frac{1}{r} = Acosh(\alpha\theta + B) + \frac{\beta}{\alpha}$, where $A$ and $B$ are determined by the initial conditions.

If $\alpha = 0$, the solution has the form $\frac{1}{r} = \beta \theta^2 + A\theta + B$

If $\alpha < 0$, $\frac{1}{r} = Acos(\alpha\theta + B) + \frac{\beta}{\alpha}$.

Answer 3

The key fact here is a very nice one that I learned from a mathematical physics book somewhere, but deserves to be better known: If $F(r)$ is a rotationally symmmetric force (with positive $F$ meaning towards the origin), and we have an almost circular orbit of radius $r_0$, then the angle between perihelions of the orbit will be $$\frac{2 \pi}{\sqrt{3+\tfrac{F'(r_0) r_0}{F(r_0)}}}.$$ When $F$ is the inverse square, we have $r F'/F = -2$, so perihelions are precisely $2 \pi$ apart, in other words, orbits close up. To solve your problem, you can look up how fast the orbit of Mercury precesses and match that with the proposed inverse cube correction term.

Derivation Let's suppose that our orbit has angular momentum $L$. That means its angular velocity is $L/(m r^2)$, so the centripetal force is $m r \left( \tfrac{L}{m r^2} \right)^2 = \tfrac{L^2}{m r^3}$. In other words, we have $$m r'' = -F(r) + \frac{L^2}{m r^3}.$$ If we had a perfect circular orbit at radius $r_0$, we would have $F(r_0) = \tfrac{L^2}{m r_0^3}$.

Now, suppose that we perturb our perfect circular orbit a little, to $r_0+\epsilon(t)$. If $\epsilon$ is small enough that we can approximate formulas with their Taylor series up to linear order, then $$m \epsilon'' = -F(r_0) - F'(r_0) \epsilon + \frac{L^2}{m r_0^3} - \frac{3 L^2}{m r_0^4} \epsilon = - \left(F'(r_0)+ \frac{3 L^2}{m r_0^4} \right) \epsilon.$$ So, by the standard analysis of a mass on a spring, $\epsilon$ oscillates with period $$2 \pi \sqrt{\frac{m}{F'(r_0) + \tfrac{3 L^2}{m r_0^4}}}.$$

To get the angle traveled between oscillations, multiply the period by the angular velocity $L/(m r_0^2)$ to get $$\theta = 2 \pi \sqrt{\frac{m}{F'(r_0) + \tfrac{3 L^2}{m r_0^4}}} \frac{L}{m r_0^2}=\frac{2 \pi}{\sqrt{3+\tfrac{m F'(r_0) r_0^4}{L^2}}}.$$ Finally, remember that $F(r_0) = \tfrac{L^2}{m r_0^3}$ to rewrite this as $$\frac{2 \pi}{\sqrt{3+\tfrac{ F'(r_0) r_0}{F(r_0)}}}.$$

Answer 4

Question 1:

This does not change qualitatively anything in the sense that your two-body problem accepts, formally speaking, the same kind of solutions than the Kepler's two body problem, i.e. conics but now paramterized by $\alpha$.

The two-body problem is best solved in the frame of the center of mass of the system. In that reference frame it appears as central force field problem applied to a particle of reduced mass $\mu = \frac{Mm}{m+M}$. Since in our case $M\gg m$, by asumption, the motion of the mass M around the center of mass can be safely ignored (the center of mass coincides with the star). Under this approximation the reduced mass $\mu\sim m$ and the total angular momentum of the two body system $\vec{L}=\mu \vec{r}\times \dot{\vec{r}}$ with that of the planet $\vec{L}\sim m \vec{r} \times \dot{\vec{r}}$.

The central problem in the fixed star reference frame has the force field given by, $$\vec{F} = \frac{-GMm}{r^2} \hat{r} +\alpha\frac{L^2}{m^2r^3} \hat{r} = -\frac{dU}{dr} \hat r$$ deriving from potential energy, $$U(r) = \frac{GMm}{r}-\frac{1}{2}{\alpha}\frac{L^2}{m^2r^2}$$

It is a well known result that in a central force field the angular moment $L$ is a first integral of the motion, i.e. stay constant, $$\dot{\vec{L}}=0.$$ This implies that the motion at any time keeps being planar, in a plane orthogonal to the angular momentum. Making use of polar coordinates in this plane of normal unit vector $\hat{e}_z$, one has, ${L} = L \hat{e}_z = mr^2\dot \phi \hat{e}_z$, yielding, $$L=mr^2\dot \phi = \mbox{Cte}.$$

Another first integral of the motion is the total energy $E$ of the system, defined as the sum of the kinetic energy and the potential energy of the planet,

$$E=\frac{1}{2} m(\dot r^2 + r^2\dot \phi^2)+U(r)=\frac{1}{2} m {\dot r^2} + \frac{L^2}{2mr^2}+U(r)$$.

This equation of energy, shows that the motion in a central field is formally similar to a "linear or translational" motion (i.e. with a single degree of freedom, the radius $r$) but in a field of "effective" potential"energy,

$$U_{\rm eff}(r)=\frac{L^2}{2mr^2}+U(r) = (1-\alpha)\frac{L^2}{2mr^2}+\frac{GMm}{r}$$

Now setting $L_\alpha = L\sqrt{1-\alpha}$, one gets, $$U_{\rm eff}(r)=\frac{L^2}{2mr^2}+U(r) = \frac{L_\alpha^2}{2mr^2}+\frac{GMm}{r}$$.

and a total energy

$$E=\frac{1}{2} m {\dot r^2} + U_{\rm eff}(r) = \frac{1}{2} m {\dot r^2} + \frac{L_\alpha^2}{2mr^2}+ \frac{GmM}{r}$$

Both equations (effective potential or total energy) describe the motion of a classical two-body Keplerian system that would have $L_\alpha$ instead of $L_0=L$ as total angular momentum. Hence we are led to conclude that a substitution of $L$ by $L_\alpha$ in all the results of the classical Keplerian two body system theory yields the correct solution to your problem. One obtains the same kind of motion than for a classical two-body Keplerian system (the case $\alpha =0$) but the trajectories instead of being entirely governed by the pair $(E,L_0)$ are now governed by $(E, L_\alpha)$. Solutions for $\alpha\neq 0$ are deduced from the Keplerian case $\alpha=0$ homothetically by the ratio $\sqrt{1-\alpha}$ (possibly pure imaginary complex when $\alpha>1$).

With a gravitational force field (of potential $\frac{GMm}{r}$) ---and more generally speaking in every Coulombian force field (of potential $-k/r$) one can easily integrate the trajectory as a conics of equation, $$p/r=1+e\cos\phi$$

where,

$$p=\frac{L^2}{km}=-\frac{L^2}{GMm^2}$$ is known as the parameter and $e$ the excentiricty of the conics is given by, $$e=\sqrt{1+\frac{2EL^2}{mk^2}}=\sqrt{1+\frac{2EL^2}{G^2M^2m^3}}$$.

Substitute $L \rightarrow L_\alpha$ in this classical result and you demonstrate that in your case, you obtain a (family of) conics (indexed by $\alpha$) of corresponding parameter(s) $p_\alpha$ and excentricitie(s) $e_\alpha$.

When $\alpha\leq 1$ the centrifugal term $\frac{L_\alpha^2}{2mr^2}$ in the effective potential is positive. (i) if $E<0$, the motion is bounded, $e_\alpha$ < 1 and you get an ellipse; (ii) if $E=0$ you get a round circle. (iii) if $E>0$, the motion is unbounded, $e_\alpha$ > 1 and you get an hyperbola. These results still hold for $\alpha=1$, simply it is the Keplerian case of null angular momentum giving a degenerated ellipse, i.e. a motion accelerated along a straight line towards the central star.

When $\alpha>1$ the centrifugal term $\frac{L_\alpha^2}{2mr^2}$ in the effective potential becomes negative. I do not see any substantial changes to the preceding results except that now, one get bounded motion along ellipses when $E>0$ and unbounded motion along hyperbolas when $E<0$.

Question 2:

I bet that $\alpha=0$ since the inclusion of your extra force term is not really necessary. Suppose you introduce a force term with let say some $L$ and $\alpha=0.1$. Your model is equivalent in a sense to one without extra force term but with planet momentum $L$ decreased to an adequate value i.e. multiplied by $\sqrt{1-\alpha}$. As a physicist, I would prefer this simpler formal solution.