Why is the temperature zero in the ground state?

Mad props for a cool question. I'm going to justify essentially the converse of the statement because it doesn't make much sense to talk about the temperature of a system that is in a pure state.

Let's assume that we're talking about a quantum system with discrete energy spectrum (with no accumulation points) in thermal equilibrium. Let $\beta = 1/kT$ be the inverse temperature. Then recall that the Boltzmann distribution tells us that the population fraction of systems in the ensemble corresponding to energy $E_i$ is given by $$ p_i = \frac{g_ie^{-\beta E_i}}{Z} $$ where $g_i$ is the degeneracy of the energy level. In particular, note that the relative frequency with which energies $E_i$ and $E_j$ will be found in the ensemble is $$ p_{ij}(\beta) = \frac{g_ie^{-\beta E_i}}{g_je^{-\beta E_j}} = \frac{g_i}{g_j}e^{-\beta(E_i - E_j)} $$ In particular, let $i=0$ correspond to the ground level, then the frequency of any level relative to the ground level is $$ p_{i0}(\beta) = \frac{g_i}{g_0}e^{-\beta(E_i - E_0)} $$ Notice that since the ground level has the lowest energy by definition, we have $E_i - E_0 \geq 0$, but zero temperature corresponds to the limit $\beta \to \infty$, and we have $$ \lim_{\beta \to \infty } p_{i0}(\beta) = \delta_{i0} $$ In other words, at zero temperature, every member of the ensemble must be in the ground energy level; the probability that a system in the ensemble will have any other energy becomes vanishingly small compared to the probability that a member of the ensemble has the lowest energy.

Addendum - September 18, 2017

In response to a question in the comments about whether or not at zero temperature the system is in a pure state:

Recall that a quantum state (density matrix) $\rho$ is said to be pure if and only if $\rho^2 = \rho$. We now show that as $T\to 0$, or equivalently as $\beta\to+\infty$, the thermal density matrix $\rho$ approaches a density matrix $\rho_*$ that is pure if the ground level is non-degenerate and not pure otherwise. We will rely on an argument quite similar in character to the one given above in which we compared the probabilities of finding a system in a given energy level when we approach zero temperature.

For any positive integer $d$, let $I_d$ denote the $d\times d$ identity matrix. As above, we consider a system with discrete energy levels $E_0<E_1<\dots$ and with corresponding degeneracies $g_0, g_1, \dots$. Recall that the thermal density matrix, namely the density matrix for a system in thermal equilibrium with a heat bath at inverse temperature $\beta$, is given in the energy eigenbasis by:

\begin{align} \rho = \frac{1}{Z} \begin{pmatrix} e^{-\beta E_0}I_{g_0} & & & \\ & e^{-\beta E_1} I_{g_1} & & \\ && e^{-\beta E_2} I_{g_2} \\ & & &\ddots \end{pmatrix}, \qquad Z = \sum_j g_j e^{-\beta E_j}. \end{align}

Let $i\geq 0$ be given, and let us concentrate on the scalar factor in front of the identity matrix $I_{g_i}$ in the $i^\mathrm{th}$ block of the density matrix:

\begin{align} \frac{e^{-\beta E_i}}{Z} = \frac{e^{-\beta E_i}}{\sum_j g_je^{-\beta E_j}} = \frac{e^{-\beta(E_i -E_0)}}{\sum_jg_j e^{-\beta(E_j - E_0)}} = \frac{e^{-\beta(E_i -E_0)}}{g_0 + \sum_{j>0}g_j e^{-\beta(E_j - E_0)}}. \end{align}

Therefore, we have

\begin{align} \lim_{\beta\to +\infty} \frac{e^{-\beta E_i}}{Z} = \frac{\lim_{\beta\to +\infty}e^{-\beta(E_i -E_0)}}{g_0 + \lim_{\beta\to +\infty}\sum_{j>0}g_j e^{-\beta(E_j - E_0)}} = \frac{\delta_{i0}}{g_0 + 0} = \frac{\delta_{i0}}{g_0}. \end{align} Care may need to be taken in the case of an infinite-dimensional Hilbert space in asserting that the limit of the sum in the denominator approaches zero since it's an infinite series. It follows that \begin{align} \rho_* = \lim_{\beta\to+\infty} \rho = \frac{1}{g_0}\begin{pmatrix} I_{g_0} & & & \\ & 0 & & \\ && 0 \\ & & &\ddots \end{pmatrix} \end{align}

and therefore in particular $\rho_*^2 = \rho_*$ if and only if $g_0 = 1$.