How to compute classical probability distribution for 1D harmonic oscillator with $K/x$ (central force) potential energy?

Question

I am trying to find, or derive, the probability distribution function for a classical 1D harmonic oscillator with a $K/x$ potential energy (from a $K/x^2$ central force). I am familiar with the derivation of the classical probability distribution for the 1D oscillator with a $kx^2/2$ potential energy, where the probability is $$P=\frac1{\pi\sqrt{x_0^2 - x^2}}$$ and +/- $x_0$ are the turning points.

Is there already an existing, well-known classical probability function solution for this $K/x$ variation on the oscillator?

I'm guessing the solution is something like $$P \approx \frac1{\sqrt{K/E-1/x}}$$ but my attempts so far to derive it (using the same approach as for the $kx^2/2$ oscillator) have not worked – doing the integral $$\int_0^{K/E}\frac1{\sqrt{E-K/x}}\;\mathrm dx$$ to determine the oscillator period does not seem to work as I'd hoped.

A motivation for this question is, for example, to formulate a linear oscillator approximation to a highly eccentric elliptical orbit about a central attracting force (e.g. gravitating mass, not necessarily on the microscopic scale). Although one can anticipate many problems arising as $r->0$ (relativistic and quantum mechanical), it seems they'd wash out on both sides of $r=0$ to yield some well characterized periodic, oscillating motion with a period defined in terms of K and m, and a classical probability distribution in some way proportional to 1/|velocity|. Is this reasonable to expect?

Thanks for suggestion on how to derive, and/or an existing solution.

Answer 1

So it sounds like you are trying to compute a classical trajectory’s “staying time” on any particular interval, as its probability of being observed in that interval. This staying time is proportional to the inverse of velocity, so the energy equation says $$ U_\text{max} = U(x) + \frac12 m v^2\\ v =\sqrt{\frac2m \big(U_\text{max} - U(x)\big)}.$$ If I am understanding your project correctly, then the only thing that you've got wrong is a sign error. Your potential is $-k/|x|$ and so $U_\text{max}$ is negative for bound states, so your probability density function looks instead like $$ \rho(x)\propto \frac{1}{\sqrt{1/|x| - 1/|x_0|}}.$$This has some pathologies but it doesn't look too bad. In particular we can multiply both top and bottom by $\sqrt{|x_0|}$ to find $$ \rho(x)\propto \frac{1}{\sqrt{|x/x_0|^{-1} - 1}},$$ suggesting that after a change of variables in the integral there is a universal cumulative density function $F(x)$ such that $F(0)=\frac12, F(1) = 1, F(-x) = 1 - F(x),$ as our boundary conditions, and because it is a CDF,$$\rho(x)=\frac{ \mathrm d}{ \mathrm dx} ~F(x/|x_0|).$$ So in some sense we can kind of just set $x_0=1$ and substitute it back in later.

What remains is to find the prefactor, so assume $x$ positive and then $$\begin{align} F(x) &=\frac12 + \alpha \int_0^x \mathrm dx' \frac{1}{\sqrt{1/x' - 1}}\\ &= \frac12 + \alpha \left[-x'\sqrt{\frac1{x'}-1}-\tan^{-1}\sqrt{\frac1{x'}-1} \right]_0^x \\ &= \frac12 + \alpha \left[\frac\pi2-x\sqrt{\frac1x-1}-\tan^{-1}\sqrt{\frac1x-1} \right] \end{align} $$ So I think $\alpha$ is just $1/\pi$? And so the full function would be $$\rho(x)=\frac1{\pi |x_0|}\frac1{\sqrt{|x/x_0|^{-1} -1}},$$or so.

Answer 2

To find the probability density you need to have an intuitive understanding of the motion under the force. For the harmonic oscillator, the particle is always oscillating from $x=-A$ to $x=+A$. Each cycle is identical to the previous one, and so the probability of finding the particle between $x$ and $x+dx$ is $dt/T$ where dt is the time the particle takes to move from $x$ to $x+dx$ and $T$ is the total time period of one oscillation. Since the particle moves the slowest near the amplitudes, it will spend a larger fraction of the total period there compared to other points, and so you would expect the probability density to be highest at those points. Which is the case, as you can see from the equation you wrote. So you were right on that point - the classical probability density should be proportional to $1/v$.

Now, in your case the potential is $k/x$, which is an attractive potential towards the center. The force gets infinitely large as you approach the center. So there is no oscillation here. So your anticipation that there is "some well characterized periodic, oscillating motion with a period defined in terms of $k$ and $m$" was wrong. The particle starts with some initial velocity $v_0$ at some initial position $x_0$. If the initial velocity is towards the center, then it approaches the center and then stops there. If the initial velocity is away from the center, then it moves away at first and then turns and approaches the center and stops there.

Nevertheless, it is possible to find the probability density if you properly define the problem. You would expect the probability density to be a Dirac delta function because the particle spends an infinite amount of time at the center and a negligible amount of time (compared to infinity) at other points. However, if you want the probability density of finding the particle during the course of it's trajectory from $x_0$ to $x=0$, then first find the time it takes to go there as a function of $x$, $$T = \int dt = \int_{x_0}^0 \frac{1}{v}dx$$ $$ = \lim_{x \to 0} \int_{x_0}^x \pm \frac{1}{\sqrt{(2/m)(E - V(x))}}$$ where I used $E = mv^2/2 + V(x)$. The value of $E$ is given by the initial conditions $$E = \frac{1}{2}mv_0^2 + V(x_0)$$ The limit sign is there because I don't know what the integral comes out to be, but since the potential goes to infinity as $x \to 0$, just putting $0$ in the upper limit might give absurd results. The plus minus sign is related to the initial condition. If the initial velocity is to the right, then you choose the plus sign. If on top of that the initial position is at positive $x$, then you break the integral into two parts. You choose the plus sign for the integral from $x_0$ to $x^*$, the turning point. And then you choose the minus sign for the integral from $x^*$ to $0$, and add up both the integrals. Similar considerations apply when the initial velocity is to the left.

Now after finding $T$, the rest is easy. $\rho dx$, which is the probability of finding the particle between $x$ and $x+dx$ is also $dt/T$. And so $$\rho dx = \frac{dt}{T}$$ $$\rho = \frac{1}{T} \frac{1}{|dx/dt|} = \frac{1}{|v|T}$$ $$\rho(x) = \frac{1}{T \sqrt{(2/m)(E - V(x))}}$$