In chapter 6 of "Charge and Energy Transfer Dynamics in Molecular Systems", the authors consider diatomic molecules in contact with a heat bath. We have established before that the nuclear coordinates $R$ are quantum mechanical and can be described in the adiabatic (Born-Oppenheimer) approximation. In this approximation, the nuclear motion is governed by the Hamiltonian
$$ H_{nuc} = T_{nuc}(R)+V_{nuc-nuc}(R)+E_a(R),$$
where $E_a$ is the ($R$-dependent) eigenvalue of the electronic Hamiltonian,
$$ \big(T_{el}+V_{el-el}+V_{nuc-el}(R) \big)\phi_a(r,R)=E_a(R)\phi_a(r,R).$$
Usually, the potential term in $H_{nuc}$ has a local minimum and can be Taylor expanded to yield a harmonic oscillator problem. The nuclear eigenfunctions can be approximated by the eigenfunctions of the harmonic oscillator.
Now the authors claim that if the temperature is much larger than the oscillator energy $k_B T >>\hbar \omega$, the molecule coordinates can be treated classically, which means in this case that they are localized at a certain value $R$, following the statistical distribution
$$P(R)=\frac{1}{Z} e^{-(V_{nuc-nuc}(R)+E_a(R))/k_B T} \tag{1}$$
(This also happens in a lot of molecular dynamics simulations: The molecular coordinates are localized and follow classical trajectories.)
Question: Why are the nuclear coordinates localized like this?
Naively, I would have expected that the nuclear system is described by the reduced density operator (Gibbs ensemble) $\rho=\frac{1}{Z} e^{-H_{nuc}/k_B T}$. There, the nuclei would occupy eigenstates of $H_{nuc}$, with the probability depending on the energy of that state. Clearly the eigenstates are delocalized across the whole harmonic potential, especially for high energy states.
Edit: I would like to point out another disparity between the QM and classical descriptions: The statistical operator $\rho$ does not evolve in time since $[H_{nuc},\rho]=0$ (neglecting non-adiabatic processes). The classical coordinates however oscillate with frequency $\omega$.