It appears that many people have tried to show adiabatic theorem from Liouville's theorem, e.g., Li's note, or at least tried to find some relations, e.g., Rugh, Adib and Tong's lecture notes Sec. 4.6.1 page 116. Can this be done rigorously? I am particularly interested in the classical (possibly statistical) mechanics sense. Thank you in advance!

Below are some random observations that may or may not be helpful in addressing the question.

An interesting observation is that the bulk entropy in statistical mechanics is an adiabatic invariant. In quantum mechanics, the integral of adiabatic invariant roughly gives the quantum number and the Planck constant, $2 \, \pi \, \hbar$. \begin{align} 2 \, \pi \, J &= \oint p \, dq \\ &= (n + \delta) \, 2 \, \pi \, \hbar, \end{align} where, $\delta$ is around $1/2$. For an almost closed system, since the quantum number $n$ can be written as a function of energy $E$ as $n(E)$, the adiabatic invariant $J$ is also a function of energy $E$.

This means that the adiabatic invariant $J(E)/\hbar$ gives the number of states $\Omega(E)$ with energy no greater than $E$, and this number is a constant. Generally, it is shown by Hertz (article 1 and article 2) that the phase-space volume enclosed by the energy surface is an adiabatic constant, a fact appreciated by Einstein. In this demonstration, Hertz used Liouville's theorem.

Now in statistical mechanics, $\Omega(E)$ is the bulk density of states[1], and it is related to the usual surface density of states $\omega(E) = e^{ S(E)/k_B }$ as \begin{align} \omega(E) = \frac{ \partial \Omega(E) } { \partial E } \Delta E. \end{align} So we may say the bulk density of states, or the bulk entropy, is an adiabatic invariant. This is, I suppose, why in thermodynamics, we say that in a reversible adiabatic process, the entropy change is zero, $dS = 0$ (or no heat is absorbed or released $dQ = 0$). According to Hertz and Rugh, we are referring to the bulk entropy.

Hertz and Rugh's proof of adiabatic invariance of the bulk entropy

The classical version is given in Hertz, the modern version is presented in Rugh.

The basic argument is the following. Recall from definition, $\Omega(E) = \int \Theta(E - H(p,q, \lambda)) \, dp \, dq$, and $\omega(E) = \Delta E \, \int \delta(E - H(p,q, \lambda)) \, dp \, dq$. Then \begin{align} d \Omega(E, \lambda) &= (\partial \Omega/\partial E) \, dE + (\partial \Omega/\partial \lambda) d\lambda \\ &= \omega(E) \, dE - \omega(E) \langle \partial H/\partial \lambda \rangle_E d\lambda, \end{align}

If $\lambda$ is changed slowly, then the time average $\overline{ dE/d\lambda } = \overline{ \partial H/ \partial \lambda }$ should match the phase space average $\langle \partial H/\partial \lambda \rangle_E$ over the isoenergetic shell, and the right-hand side vanishes.

References

Hans Henrik Rugh, A Micro-Thermodynamic Formalism, nlin/0201062.

P. Hertz, Ann Phys. (Leipzig) 33, 537 (1910). Particularly pages 546 to 548. Notations differences: energy: $E \rightarrow \varepsilon$; parameter: $\lambda \rightarrow a$; Bulk density of states: $\Omega \rightarrow V$; surface density of states, $\omega \rightarrow \omega$; breadth of the energy shell, $\Delta E \rightarrow d\nu \sim \Delta E/|\nabla E|$.
[1] Some people define $k_B \log \Omega(E)$ as the “bulk entropy”, and this leads to some interesting debate, such as Dunkel and Hilbert and Frenkel and Warren.