Bertrand's theorem and nearly-circular motion in a Yukawa potential

In the other post, I explained why the two solutions are both wrong. Here, I try to derive the correct answer of: $$ piBig(frac{r_0}{a}Big)^2 $$ To arrive at the actually correct answer, we first find the condition for having a circular orbit. To accomplish this end, we use exactly the force law that you point out: $$ f(r) = -frac{k}{r^2}e^{-frac{r}a} - frac{k}{ra}e^{-frac{r}{a}} $$

When orbits are circular, centripetal acceleration times mass is equal to force, or: $$ begin{align*} -frac{k}{r_0^2}e^{-frac{r_0}{a}} - frac{k}{r_0a}e^{-frac{r_0}{a}} &= -mr_0dot{theta}^2 &=-frac{l^2}{mr_0^3} end{align*} $$ where $l = mr^2 dot{theta}$. An orbit is circular if and only if the condition above is fulfilled. Suppose such a solution exists for the given values of $k$ and $a$ and suppose that this orbit is perturbed very slightly to $r(t) = r_0 + epsilon(t)$ where $epsilon(t)$ is a function that represents the small difference between the radius of the circular orbit and the actual orbit at time $t$. It has a local maximum at $t = theta = 0$, where accordingly $frac{mathrm{d}}{mathrm{d}t} epsilon = 0$ and $epsilon(t=0) > 0$

Let us first attempt the obvious method, Newton's Laws. As we will see, this fares very poorly, and we will ultimately turn to the Jacobi integral (conservation of energy) instead. (So if you want to see the full solution, skip the next section, and go down to the one below it)

Newton's second law tells us that: $$ begin{align*} m cdot a_{mathrm{radial}} &= F_{mathrm{radial}} mddot{r}-mrdot{theta}^2 &= -frac{k}{r^2}e^{-frac{r}a} - frac{k}{ra}e^{-frac{r}{a}} m{epsilon''(t)} - frac{l^2}{m(r_0+epsilon(t))^3} &= -frac{k}{(r_0 + epsilon(t))^2}e^{-frac{r_0+epsilon(t)}a} - frac{k}{(r_0 + epsilon(t))a}e^{-frac{r_0 + epsilon(t)}{a}} end{align*} $$ From now on, for the sake of sanity, we write $epsilon:= epsilon(t)$, and we also factor out $r_0$'s from each denominator: $$ begin{align*} mddot{epsilon} - frac{l^2}{mr_0^3(1+frac{epsilon}{r_0})^3} &= -frac{k}{r_0^2(1 + frac{epsilon}{r_0})^2}e^{-frac{r_0+epsilon}a} - frac{k}{ar_0(1 + frac{epsilon}{r_0})}e^{-frac{r_0 + epsilon}{a}} end{align*} $$ Next, we substitute the circular orbit condition above into the LHS, where we find: $$ begin{align*} mddot{epsilon} + Big(-frac{k}{r_0^2}e^{-frac{r_0}{a}} - frac{k}{r_0a}e^{-frac{r_0}{a}}Big)frac{1}{(1+frac{epsilon}{r_0})^3} &= -frac{k}{r_0^2(1 + frac{epsilon}{r_0})^2}e^{-frac{r_0+epsilon}a} - frac{k}{ar_0(1 + frac{epsilon}{r_0})}e^{-frac{r_0 + epsilon}{a}} &= -ke^{-frac{r_0}{a}}Big(frac{e^{-frac{epsilon}a}}{r_0^2(1 + frac{epsilon}{r_0})^2} + frac{e^{-frac{epsilon}a}}{ar_0(1 + frac{epsilon}{r_0})}Big) end{align*} $$ Now, expand this expression using the binomial approximation, avoiding the mistake that Reina makes; because we intend on square rooting our approximation down the line, we must expand everything to second order: $$ begin{align*} mddot{epsilon} &+ Big(-frac{k}{r_0^2}e^{-frac{r_0}{a}} - frac{k}{r_0a}e^{-frac{r_0}{a}}Big)(1 - 3 frac{epsilon}{r_0} + 6frac{epsilon^2}{r_0^2}) &= -ke^{-frac{r_0}{a}}Big(frac{e^{-frac{epsilon}a}}{r_0^2}{(1 - 2frac{epsilon}{r_0} + 3frac{epsilon^2}{r_0^2})} + frac{e^{-frac{epsilon}a}}{ar_0}{(1 - frac{epsilon}{r_0}+ frac{epsilon^2}{r_0^2})}Big) end{align*} $$ And expand $e^{-frac{epsilon}{a}} approx 1 - frac{epsilon}{a} + frac{1}{2} frac{epsilon^2}{a^2}$ and substitute in the equation above. Then, carefully expand, and drop all terms of $mathcal{O}(epsilon^3)$ and higher. When all is said and then, the result should be: $$ begin{align*} 0 &= m ddot{epsilon} &+ frac{k}{r_0^2}e^{frac{-r_0}{a}}(frac{3 epsilon}{r_0} - frac{epsilon}{a} + frac{1}{2}frac{epsilon^2}{a^2} + 2 frac{epsilon^2}{r_0a} - 3frac{epsilon^2}{r_0^2}) &+ frac{k}{ar_0}e^frac{-r_0}{a}(frac{3epsilon}{r_0} - frac{epsilon}{a} + frac{1}{2}frac{epsilon^2}{a^2} + frac{epsilon^2}{r_0a} - 5frac{epsilon^2}{r_0^2}) end{align*} $$ which is a second-order, non-linear ordinary differential equation whose full solution can only be expressed in terms of a hypergeometric series. One can attempt to approximate this by expanding $epsilon$ into a second-order Fourier series: $epsilon(t) approx a_1 cosbeta t + a_2 cos 2beta t$ and substituting this into the equation, and equating both sides of the diff eq. This can tell you the frequency of $epsilon$'s oscillation through $beta$. I highly recommend against doing this, as it will clearly become very messy, very quickly.

Instead, when in doubt, always try conservation of energy as well. Here, this equation is: $$ E = frac{m}{2}v^2 - frac{k}{r}e^{-frac{r}{a}} $$ Recall that velocity can be split into two components, radial and orthogonal: $vec{v} = dot{r}hat{r} + rdot{theta}hat{theta}$. So, we get that $v^2 = dot{r}^2 + r^2dot{theta}^2$. Once again, expand $r = r_0 + epsilon$ and substitute both of these into the energy equation to find: $$ begin{align*} E &= frac{m}{2}dot{epsilon}^2 + frac{m}{2}(r_0 + epsilon)^2dot{theta}^2 - frac{k}{r_0 + epsilon} e^{-frac{r_0 + epsilon}{a}} &= frac{m}{2}dot{epsilon}^2 + frac{m^2(r_0 + epsilon)^4dot{theta}^2 }{2m(r_0+epsilon)^2}- frac{k}{r_0 + epsilon} e^{-frac{r_0 + epsilon}{a}} &= frac{m}{2}dot{epsilon}^2 + frac{l^2}{2m(r_0+epsilon)^2} - frac{k}{r_0 + epsilon} e^{-frac{r_0 + epsilon}{a}} &= frac{m}{2}dot{epsilon}^2 + frac{l^2}{2mr_0^2(1+frac{epsilon}{r_0})^2} - frac{ke^{frac{-r_0}{a}}}{r_0(1 + frac{epsilon}{a})} e^{-frac{epsilon}{a}} &= frac{m}{2}dot{epsilon}^2 + frac{1}{2}frac{l^2r_0}{mr_0^3(1+frac{epsilon}{r_0})^2} - frac{ke^{frac{-r_0}{a}}}{r_0(1 + frac{epsilon}{a})} e^{-frac{epsilon}{a}} end{align*} $$ Once again, substitute the circularity condition: $$ begin{align*} E &= frac{m}{2}dot{epsilon}^2 + frac{1}{2}Bigg[frac{k}{r_0^2}e^{-frac{r_0}{a}} + frac{k}{r_0a}e^{-frac{r_0}{a}}Bigg]frac{r_0}{(1+frac{epsilon}{r_0})^2} - frac{ke^{frac{-r_0}{a}}}{r_0(1 + frac{epsilon}{a})} e^{-frac{epsilon}{a}} end{align*} $$ And as before, avoiding the mistake that Reina makes, as we intend on later square-rooting our approximation, expand the binomial and Taylor series to the second order: $$ begin{align*} E &= frac{m}{2}dot{epsilon}^2 + frac{1}{2}Bigg[frac{k}{r_0^2}e^{-frac{r_0}{a}} + frac{k}{r_0a}e^{-frac{r_0}{a}}Bigg]{r_0}{(1-2frac{epsilon}{r_0}+3frac{epsilon^2}{r_0^2}) } - frac{ke^{frac{-r_0}{a}}}{r_0}(1 - frac{epsilon}{a}+ frac{epsilon^2}{a^2}) (1 - frac{epsilon}{a} + frac{epsilon^2}{2a^2}) end{align*} $$ Carefully, as in very carefully, as in (do-this-twice-on-two-different-sheets-of-paper-and-verify-your-results)-level of carefully, expand this expression and drop all terms of order $mathcal{O}(epsilon^3)$ and higher. You should then get: $$ begin{align*} E &= frac{mdot{epsilon}^2}{2} + ke^{-frac{r_0}{a}}bigg[frac{1}{2}frac{1}{r_0} - frac{epsilon}{r_0^2} + frac{3}{2}frac{epsilon^2}{r_0^3} + frac{1}{2}frac{1}{a} - frac{epsilon}{r_0a} + frac{3}{2}frac{epsilon^2}{ar_0^2}bigg] &+ ke^{-frac{r_0}{a}}bigg[-frac{1}{r_0} + frac{epsilon}{r_0^2} - frac{epsilon^2}{r_0^3} + frac{epsilon}{r_0a} - frac{epsilon^2}{r_0^2a} - frac{1}{2}frac{epsilon^2}{r_0a^2}bigg] end{align*} $$ Note that all the terms with just $epsilon$ cancel out in the expression above, and we are left with only quadratics ($epsilon^2$) and constants: $$ begin{align*} E &= frac{mdot{epsilon}^2}{2} + ke^{-frac{r_0}{a}}bigg[frac{1}{2}frac{epsilon^2}{r_0^3} + frac{1}{2}frac{epsilon^2}{ar_0^2} - frac{1}{2}frac{epsilon^2}{r_0a^2}bigg] &+ ke^{-frac{r_0}{a}}bigg[frac{1}{2}frac{1}{a} - frac{1}{2}frac{1}{r_0}bigg] end{align*} $$ Call the last term on the RHS $C$ as it is constant, and subtract it from both sides, in the process defining $E^* = E - C$: $$ begin{align*} E^* = E - C &= frac{mdot{epsilon}^2}{2} + ke^{-frac{r_0}{a}}bigg[frac{1}{r_0^3} + frac{1}{ar_0^2} - frac{1}{r_0a^2}bigg]frac{epsilon^2}{2} &= frac{mdot{epsilon}^2}{2} + frac{Kepsilon^2}{2} end{align*} $$ where $K$ is the big coefficient in front of the second term. Observe that this is the energy equation for a simple harmonic oscillator. You noted here that you had some difficulties understanding how to find the apside advance from frequency. Instead, I recommend thinking of this in terms of the period, which is given for the SHO corresponding to $epsilon$ by: $$ T_mathrm{pert} = 2pi sqrt{frac{m}{K}} = 2 pi sqrt{frac{m}{ke^{-frac{r_0}{a}}bigg[frac{1}{r_0^3} + frac{1}{ar_0^2} - frac{1}{r_0a^2}bigg]}} $$ On the other hand the period of the orbit can be found by $frac{2pi}{dot{theta}}$. We can extract $dot{theta}$ from the circularity condition: $$ begin{align*} mr_0dot{theta}^2&= frac{k}{r_0^2}e^{-frac{r_0}{a}} + frac{k}{r_0a}e^{-frac{r_0}{a}} = ke^{-frac{r_0}{a}}bigg[frac{1}{r_0^2} + frac{1}{r_0a}bigg] Longrightarrow dot{theta} &=sqrt{frac{ke^{-frac{r_0}{a}}bigg[frac{1}{r_0^3} + frac{1}{r_0^2a}bigg]}{m}} Longrightarrow T_{mathrm{orbit}} &= 2 pi sqrt{frac{m}{ke^{-frac{r_0}{a}}bigg[frac{1}{r_0^3} + frac{1}{r_0^2a}bigg]}} end{align*} $$ Just by looking at the equations, we can tell that $T_{mathrm{pert}} > T_{mathrm{orbit}}$ as the denominator in the former is less than in the latter.

Let's think about an analogy for this situation: if it takes Bob $T_B = 5$ minutes to paint a wall and Alice $T_A = 3$ minutes, and if they start at the same time, and continuously paint as many walls as they can, when the slower person (Bob) has finished his first wall, how many additional walls has Alice finished? Clearly, its $$ frac{2}{3} mathrm{walls} = frac{5 mathrm{min} - 3 mathrm{min}}{3min} = frac{T_B - T_A}{T_A} = frac{T_B}{T_A} - 1 $$ This is just true in general. Here, the perturbation is like Bob, it is slower than the orbital motion, Alice. Therefore, the number of additional (more than $1$) cycles the orbit has completed by the time the perturbation undergoes its first cycle is: $$ begin{align*} frac{T_{mathrm{pert}}}{T_{mathrm{orbit}}} -1 &= sqrt{frac{frac{1}{r_0^3} + frac{1}{r_0^2a}}{frac{1}{r_0^3} + frac{1}{ar_0^2} - frac{1}{r_0a^2}}} - 1 &= sqrt{frac{1 + frac{r_0}{a}}{1 + frac{r_0}{a} - frac{r_0^2}{a^2}}} - 1 &= sqrt{frac{1+x}{1+x - x^2}} - 1 (mathrm{where} x := frac{r_0}{a}) &=sqrt{frac{1+x-x^2}{1+x-x^2} + frac{x^2}{1+x-x^2}} - 1 &= sqrt{1 + frac{x^2}{1 + x - x^2}} - 1 &approx1 + frac{1}{2}bigg(frac{x^2}{1 + x - x^2}bigg) - 1 &=frac{1}{2}x^2bigg(frac{1}{1+(x - x^2)}bigg) &= frac{1}{2}x^2(1 - (x-x^2) + mathcal{O}((x-x^2)^2)) &= frac{1}{2}x^2 + mathcal{O}(x^3) = frac{1}{2}frac{r_0^2}{a^2} + mathcal{O}(x^3) end{align*} $$ Dropping all terms of high order $mathcal{O}(x^3)$, we find that the orbit has underwent an additional $frac{1}{2}frac{r_0^2}{a^2}$ cycles or that we are now: $$ 2 pi cdot frac{1}{2}frac{r_0^2}{a^2} = pibig(frac{r_0}{a}big)^2 $$ radians into the second cycle of the orbit. Thus, the peak distance away from the origin (which obviously occurs when $mathrm{pert}$ has undergone a full cycle) is $pi(frac{r_0}{a})^2$ of the way into the orbit, as desired.

Addendum: To address @YamanSanghavi's question. There are two necessary premises for the solution above:

1. Circular orbits
2. Stability of the orbit

The former has been found by setting centripetal force equal to mass times centripetal acceleration: $$ underbrace{-frac{l^2}{m} + kre^{-frac{r}{a}} + frac{kr^2}{a}e^{-frac{r}{a}}}_{text{both sides multiplied by }r^3} = 0 $$ The latter is equivalent to the positive-definiteness of the effective potential; if positive definite, the potential resembles an upward facing parabola near the local minimum, and so, stability is achieved. Here, $V_{mathrm{eff}}$ can be found via the Jacobi integral: $$ begin{align*} frac{mdot{r}^2}{2} &+ {frac{mr^2dot{theta}^2}{2} - frac{k}{r}e^{-frac{r}{a}}} = E frac{mdot{r}^2}{2} &+ underbrace{frac{l^2}{2mr^2} - frac{k}{r}e^{-frac{r}{a}}}_{text{effective potential}} = E end{align*} $$ Its positive-definiteness can be examined second derivative by taking the second derivative with respect to $r$: $$ frac{mathrm{d}^2V_{mathrm{eff}}}{mathrm{d}r^2} = frac{3l^2}{mr^4} - frac{2ke^{-frac{r}{a}}}{r^3} - frac{2ke^{-frac{r}{a}}}{ar^2} - frac{ke^{-frac{r}{a}}}{a^2r} $$ Substitute $l^2$ from the earlier circularity condition into the second derivative. Expanding the expression for the second derivative should yield: $$ {r^3e^{frac{r}{a}}}frac{mathrm{d}^2V_{mathrm{eff}}}{mathrm{d}r^2} = k + kfrac{r}{a} - kfrac{r^2}{a^2} $$ Treating the RHS as a function of $frac{r}{a}$, we see it is the expression for a downward facing parabola. The sum of its roots is positive but their product is negative, via Vieta's formulae. Try drawing a graph of such an expression to see that the parabola's output is positive only in the case that $frac{r}{a}$ is less than the function's greater root. We find this root via the quadratic formula:

$$ frac{r}{a} < frac{-1 + sqrt{5}}{2} $$ which is a sufficiently tight upper bound on $x$, for all practical purposes.