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

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.

• Hello! I have edited your question using MathJax (LaTeX) math typesetting. For future questions, you can refer to MathJax basic tutorial and quick reference. Thanks! Aug 4, 2021 at 18:15
• If you change the potential from $0.5kx^2$ to $k/x$ you no longer have a harmonic oscillator. The $1/x$ term is closer to the radial equation for the hydrogen atom. You could see Griffiths's Intro to Quantum Mechanics for that solution. If you set the angular momentum quantum number $l$ to zero and relabel variables you should get your desired solution. Aug 4, 2021 at 18:22
• I didn't realize that the term "harmonic oscillator" applied only to a 0.5kx^2 potential, but it seems that for the k/x potential there would still be periodic linear (zero angular momentum) oscillation about the origin, correct? I do not have access to Griffith's book, but does the approach you suggest in fact yield a probability curve that is roughly P=~1/sqrt(K/E - 1/x) ? Sorry I'm not able to deduce it from the information provided.
– Dave
Aug 4, 2021 at 19:49
• "but it seems that for the k/x potential there would still be periodic linear (zero angular momentum) oscillation about the origin, correct?" No. You're looking at something like the 1D hydrogen atom: royalsocietypublishing.org/doi/10.1098/rspa.2015.0534
– Gert
Aug 4, 2021 at 20:31
• Thanks for that reference. I can see how that would be similar, but I'm actually just thinking of this in classical (non-microscopic) terms, for example as a linear oscillator approximation to a highly eccentric elliptical orbit about a central gravitating mass. Although one can anticipate many problems arising as r->0 (e.g. 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|.
– Dave
Aug 4, 2021 at 22:16

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.

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))}}$$

• I'm working through your answer and do not understand the statement "The force gets infinitely large as you approach the center. So there is no oscillation here." I understand there are several relativistic (and quantum) effects to ignore in this problem, but as the particle approaches x=0, doesn't its momentum also become "just infinite enough" to overcome the infinite attractive force and emerge on the other side with a steadily diminishing velocity that is zero at the other turning point (consider something like the limit of an orbital ellipse as eccentricity approaches zero)? Thx.
– Dave
Aug 9, 2021 at 11:15
• If the particle has no angular momentum, then it will just collide with the gravitiating mass. If there is some angular momentum, it will induce a term $\sim L/r^2$ which will stabilize the potential near the origin and create the behavior of an orbit, or at least a slingshot. But your potential just describes a particle heading straight radially toward a gravitating mass. Aug 13, 2021 at 6:55