Exact energy distribution function in time-dependent harmonic oscillator

Following a recent work by Robnik and Romanovski (J.Phys.A: Math.Gen. {\bf 39} (2006) L35, Open Syst. & Infor. Dyn. {\bf 13} (2006) 197-222) we derive the explicit formula for the universal distribution function of the final energies in a time-dependent 1D harmonic oscillator, whose functional form does not depend on the details of the frequency $\omega (t)$, and is closely related to the conservation of the adiabatic invariant. The normalized distribution function is $P(x) = \pi^{-1} (2\mu^2 - x^2)^{-{1/2}}$, where $x=E_1- \bar{E_1}$, $E_1$ is the final energy, $\bar{E_1}$ is its average value, and $\mu^2$ is the variance of $E_1$. $\bar{E_1}$ and $\mu^2$ can be calculated exactly using the WKB approach to all orders.