Boltzmann–Gibbs-distribution as resulting from a limiting density of states? I'm interested in the relation between the probability distribution $p_i$ over states of a system on the one side and the density of states $\rho(\eta)$ of its environment. (Meaning, $\int_{\eta_a}^{\eta_b} \rho(\eta) ~ \mathrm{d} \eta$ is the number of environment states with energies in the interval $[\eta_a, \eta_b]$.)
If the whole (system + environment) is energetically closed ("isolated") with a total energy $E = e + \eta$, but system and environment are in thermal equilibrium (i.e. the whole is described by the microcanonical ensemble), then it holds
$$
p_i = \frac{ \rho(E - e_i) }{ \sum_i \rho(E - e_i) }.
$$
This means, the probability distribution over states of the system is determined by
a) something that only characterizes the energetic structure of the system, the $e_i$s
b) something that only characterizes the energetic structure of the environment, $\rho(\eta)$, and
c) the total energy $E$.
This relation holds generally, for arbitrarily small or large systems and/or environments. Please note that we have not yet taken any limits!
If we now consider the thermodynamic limit, i.e. an environment composed of an infinite number of subsystems, the probability distribution $p_i$ over states of the system becomes the Boltzmann–Gibbs-distribution (aka canonical ensemble)
$$
p_i = \frac{ \exp(- \beta ~ e_i) }{ \sum_i \exp(- \beta ~ e_i) },
$$
where the sum in the denominator is called the partition function.
Using the first relation above, this distribution could now be interpreted as corresponding to a limiting density of states of the environment of the form
$$
\rho(\eta) \propto \exp( \beta \eta )
$$
which characterizes the "infinite environment". However, the expression refers to the parameter $\beta$ of the Boltzmann–Gibbs distribution, which represents the temperature and depends on the total energy $E$ (per subsystem). Whereas in the finite case $E$ only serves to connect $\rho(\eta)$ and $p_i$, it here defines $\rho(\eta)$ itself.
To me this suggests that it does not make sense to characterize an infinite environment by a density of states — but maybe there's some way around this? Or is there a mistake in the derivation somewhere else?
 A: When one proves that a small part of a greater system is described canonical ensemble, even though the greater system is described by a microcanonical ensemble, the key point is that the density of states of the greater system has the exponential form you mention, over a certain interval of energy.
Specifically what is important is that $\log \rho(\eta) = {\rm const} + \beta \eta$ over the range $\eta = \langle\eta \rangle \pm \Delta e$, where $\Delta e$ is the energy fluctuation of the small system, and $\langle \eta \rangle = E - \langle e \rangle$ is the expected energy. The value of $\rho(\eta)$ is not really important for energies that are very far outside this range, since those energies never occur.
In practice the form of $\rho(\eta)$ is usually something like $\log\rho(\eta) \approx {\rm const} + N \log \eta$ for some very large $N$ (such a form is found for gases, in particular). This log can be taylor expanded to yield the needed form for the canonical ensemble. For a large environment (one that is essentially in the thermodynamic limit and much larger than the attached system), the first order expansion is very accurate over the relevant energy range.
A: The right relationship between the two (the two distributions are pretty much the microcanonical and canonical ensemble) is that the first distribution – all states with some fixed energy – is almost exactly the same (the difference goes to zero if the environment becomes infinite) as the Boltzmann-Gibbs distribution for an appropriate inverse temperature $\beta$. 
However, the right value of $\beta$ has to be carefully adjusted for the two distributions to match in the limit. This appropriate temperature may be essentially extracted from the energy of the environment (or energy of everything, it isn't any different in the infinite-environment limit) divided by its number of degrees of freedom, $\beta=\beta(E/N)$.
A: I'm revisiting my own question after almost a year to add an answer that to my mind clears the matter up. This answer does not add much beyond @Nanite's, but spells out the details. Comments are welcome!
If the environment has an $N$-dimensional state space, for large $N$ its density of states over energy $\eta$ will typically* have the form
$$
\rho(\eta) = a_N ~ \eta ^ {N-1}
$$
($-1$ because the $\rho$ is the derivative of the state space volume).
For the numerator of the probability distribution $p_i$ it follows
$$
\rho(E - e_i) = a_N ~ (E - e_i)^{N-1} = a_N E ~ \left (1 - \frac{e_i}E \right )^{N-1}.
$$
If for increasing $N$ we keep the total energy per dimension constant, $\frac EN = \epsilon$,
$$
\rho(E - e_i) = a_N N \epsilon ~ \left (1 - \frac{e_i/\epsilon}{N} \right )^{N-1},
$$
the power term becomes $\exp(-\frac{e_i}\epsilon)$ in the limit $N \to \infty$.
However, the limit of the full state density $\rho(E - e_i)$ depends on the behavior of $a_N N$. A typical* behavior appears to be $a_N \propto \frac1{N!}$, and therefore the density converges to 0. An infinite environment can therefore not be meaningfully characterized by a density of states, even locally.
Proper limiting behavior for the probability distribution does arise since
$$
p_i = \frac{ \rho(E - e_i) }{ \sum_i \rho(E - e_i) }
= \frac{ \left (1 - \frac{e_i/\epsilon}{N} \right )^{N-1} }{ \sum_i \left (1 - \frac{e_i/\epsilon}{N} \right )^{N-1} }
$$
which becomes the Boltzmann–Gibbs-distribution
$$
p_i = \frac{ \exp(-\frac{e_i}\epsilon) }{ \sum_i \exp(-\frac{e_i}\epsilon) }
$$
for $N \to \infty$ — and we see that the temperature of the small system is determined by the energy per dimension of the environment, $1/\beta = \epsilon$.
What we can do is use this limiting behavior to locally approximate $\rho(\eta)$ around $\eta = N \epsilon$ for large but finite $N$:
$$
\rho(\eta) \approx a_N \exp(1 - N) \frac N\beta ~ \exp ( \beta ~ \eta )
\quad \propto \exp( \beta ~ \eta )
$$
The dependence on $\beta$ now appears as a natural consequence of the circumstance that we approximate locally around $N \epsilon = \frac N\beta$.

*) The weak point of this answer is of course that it assumes a "typical" state density of a high-dimensional environment of the form
$$
\rho(\eta) = a_N ~ \eta ^ {N-1}
$$
or even
$$
\rho(\eta) = a_0 ~ \frac1{N!} ~ \eta ^ {N-1}
$$
with $a_0$ not or only weakly dependent on $N$.
