In a NVE ensemble, can a particle access a state with a higher energy than the constraint? I should start with saying that I know the answer is obviously no but I am trying to make sense of it mathematically.
Let's say, for example, that each particle of a system with total energy $E$ in a NVE ensemble can have energies $ε_i$, $i=0,1,...$ with $ε_j > ε_k$ for $j > k$.
From statistical mechanics I know that the probability of finding a particle in the $i$ state is proportional to an exponential $\exp\left(-\frac{ε_i}{k_BT}\right)$. So that means that it's never $0$. Practically it is $0$ but in infinite time and high temperature would there not be a single particle that reaches a state with $ε_i > E$? What prevents this situation from happening? Is this just a problem of the mathematical model?
 A: You are mixing results from the microcanonical and canonical ensemble.  Macroscopic results derived from the two are the same, but for questions about rare events (such as enormous energy fluctuations), the two ensembles are inequivalent.
A microcanonical ensemble (which you refer to by the seldom used term "$NVE$ ensemble") is one in which the gas system has fixed values for the three extrinsic thermodynamic quantities $N$ (the number of particles), $V$ (the volume), and $E$ (the total energy).  Once these three quantities are set, every microscopic state of the system that is consistent with these three conditions is equally likely.  The usual procedure is to show that there is a most likely set of values for the microscopic energy distribution of the particles, and then that very nearly all possible states have momenta close to the most likely values.
For more complicated calculations—of things like transport properties in the gas, for example—this method is not so useful, however.  The issue is that the distribution of momentum for a single particle depends on the energy available for that particle—and hence on the energies of all the other particles, since all the particles energies have to add up to a fixed $E$.  So instead it is easier to work in the canonical ensemble, in which the total energy is allowed to fluctuate, with the temperature $T$ fixed instead of $E$.  In the canonical ensemble (and provided that the particles have negligible energies of interaction, except during their instantaneous collisions), the probability of a single-particle state $j$ with energy $\varepsilon_{j}$ is exactly $\exp\left(-\frac{\varepsilon_{j}}{kT}\right)$, and the energies of the individuals particles are independent.
When $\varepsilon_{j}$ is a reasonably likely energy, meaning ${\cal O}(kT)$ or less, the microcanonical ensemble gives essentially the same distribution.  However, when the energy is extremely large, the two ensembles give different results.  In the microcanonical ensemble, the probability of finding a particle with energy $\epsilon>E$ is simply zero, while it is $\exp\left(-\frac{\varepsilon}{kT}\right)\lesssim \exp(-N)$—which is merely incredible small—in the canonical ensemble.  For a real physical system, which answer is correct depends on what precisely you know about the system.  If you genuinely know the exact energy, then the microcanonical result is correct; if instead the temperature is known exactly, then the canonical result is correct.  Most likely, what you really know about the system is not either of those but some other closely related quantity.  And for most thermodynamic quantities, it makes no difference; whatever observable you use to characterize the system's energy, the system is overwhelmingly likely to have an essentially Maxwell-Boltzmann distribution of energies.
A: I'll add a more direct argument to Buzz's answer which already contains all the key observations.
Let's consider the case of an ideal classical gas in one dimension, for simplicity. While the canonical ensemble predicts a one particle distribution $f(\epsilon)$  proportional to $exp\left(-\frac{\epsilon}{k_BT}\right)$, where the one particle energy is $\epsilon=\frac{p^2}{2m}$, the microcanonical one-particle distribution differs for two main ingredients:


*

*it depends on the total energy $E$ as a parameter, and not on $T$;

*it must be exactly zero for $\epsilon>E$.


The actual calculation is not too complicate ( see here for instance ) and the exact result is 
$$
\begin{eqnarray}
f(\epsilon) &=& C \left(  1 - \frac{\epsilon}{E}  \right)^{\frac{(N-1)}{2}} ~~~~~~~~~~~~~~~&{\text{for}}~~~~0 \le \epsilon \le E\\
& & 0 &{\text{for}}~~~~\epsilon>E
\end{eqnarray}
$$
where $C$ is the normalization constant (depending on $E$).
Such a distribution goes to the Boltzmann's distribution at the thermodynamic limit ($E \rightarrow \infty$, $N \rightarrow \infty$, while keeping $\frac{E}{N}$ constant).
