How is quantum mechanics consistent with statistical mechanics?

The question doesn't make sense, as not every system has to be in thermal equilibrium.

The argument you're making works for classical mechanics too. "Consider a ball lying on the ground." "That's impossible, since by statistical mechanics, the position of the ball must obey the Boltzmann distribution."

You're comparing the wrong things - you must take the appropriate limits when using statistical physics (or statistical mechanics as you called it) to recover quantum and classical results. Also, recall that when we use statistical physics, we must consider the system in an ensemble formalism, i.e. the canonical ensemble. That is, we consider a (statistically) large number of identically prepared copies of our system: a single harmonic oscillator. Additionally, we must use the quantum version of statistical physics, where the ensemble is characterized by the density operator. For a nice intro to this, see Sakurai.

This article should answer your questions regarding the case of a single 1D harmonic oscillator potential viewed statistically. For a given temperature, we use the quantized energy of the oscillator to find the partition function in the canonical ensemble.

Consider the two limits:

1) Classcial: the thermal energy is much greater than the spacing of oscillator energies. Here, one recovers the classical result that the energy is $kT$.

2) Quantum: the thermal energy is much less than the spacing of oscillator energies. Here, one recovers the quantum result that for a given frequency the energy is $\frac{1}{2}\hbar \omega$

Further, using the appropriate density operator, you can recover the probabilities that you sought:

the density operator represented as a matrix in the basis $\{\lvert \phi_1 \rangle, \lvert \phi_2 \rangle \}$, for the system given by pure state $\Psi$ is

$$ \rho = \lvert \Psi \rangle \langle \Psi \lvert = \begin{pmatrix} \frac{1}{2} & \frac{1}{2} \\\ \frac{1}{2} & \frac{1}{2} \end{pmatrix} = \frac{1}{2} \begin{pmatrix}1 & 1 \\\ 1 & 1\end{pmatrix}$$

To find the probability of being in state $i$, take,

$$\mathcal{P}_i = Tr(\rho P_i)$$

where $P_i$ is the projection operator of the $i^{th}$ state, and Tr means trace.

This question does make sense.

It is true that statistical mechanics (or better most of it) is usually limited to ensembles for equilibrium systems and its main task is to provide the basis for the thermodynamic behavior of systems with many degrees of freedom. However this does not imply that it is always meaningless to speak about probabilities of states of even a single particle, provided one is interested to the average behavior in the ensemble, which will correspond to time averages, if the system made by the individual particle and the many degrees of freedom of the canonical ensemble bath correspond to a globally ergodic system.

Therefore, a single particle in interaction with a thermal bath is a perfectly admissible system on which Statistical Mechanics allows to make statements about probabilities.

The only difference, which seems to puzzle the OP, between a quantum oscillator whose state is described by the superposition of two degenerate states of energy $E_i$ and the case of the same quantum oscillator in equilibrium with a thermal bath (canonical ensemble), is that the two situations correspond to different events characterized by different external conditions. So, no surprise if the corresponding probabilities are different.

• In the case of the system in the state $\Psi$, there are only two mutually exclusive events (the system is in state $\phi_1$, or the system is in state $\phi_2$), QM allows to get the probability of the two cases with the proposed state $\Psi$, and for each case it turns out to be equal to $0.5$.
• In the case of the quantum harmonic oscillator (by the way it must be at least a 2D oscillator, if there are degenerate states) the condition of thermal equilibrium with the bath implies that not only the two degenerate states will be visited by the oscillator, but in principle the quantum oscillator will visit all the states. Each state appearing with a frequency proportional to the Boltzmann factor of the corresponding energy (the proportionality factor being the normalization of the probability distribution which is nothing but the quantum canonical partition function of the harmonic oscillator). So, it is not a surprise if the probability of the two states $\phi_1$ and $\phi_2$ are different from $0.5$.

It is interesting to notice that there is something which remains the same in the two situations: if the only degenerate states at energy $E_i$ are $\phi_1$ and $\phi_2$, even if their probability is not $0.5$ anymore, the two probabilities remain equal each other, the Boltzmann factor being the same.