Derivation of Boltzmann distribution - two questions I understand that the usual way of deriving the Boltzmann distribution involves considering a small system of energy $\epsilon$ embedded in a much larger heat bath of energy $E - \epsilon$ and the total system energy is $E$.
Given that the bath has accessible microstates $\Omega_b(E - \epsilon)$ and the system has microstates $\Omega_s(\epsilon)$ and the total system + bath has microstates $\Omega_t(E)$. 
We now state that the probability of finding the system in energy $\epsilon$ is simply 
\begin{equation}
P(\epsilon) = \frac{\Omega_s(\epsilon)\Omega_b(E - \epsilon)}{\Omega_t(E)}
\end{equation}
$\textbf{Question 1}$: This is usually replaced by 
\begin{equation}
P(\epsilon) = \frac{\Omega_b(E - \epsilon)}{\Omega_t(E)}
\end{equation}
Why should this be so? Why can we ignore the fact that the number of accessible microstates in the system that should vary with $\epsilon$?
Moving on, the usual trick is to use logarithms, since the $\Omega$ terms are all very large. Thus,
\begin{equation}
\ln P(\epsilon) = \ln\Omega_b(E - \epsilon) -\ln \Omega_t(E)
\end{equation} 
A Taylor expansion of the first term about $E$ gives us
\begin{equation}
\ln P(\epsilon) = c - \epsilon\frac{\partial\ln\Omega_b(E)}{\partial E} + O(\epsilon^2),
\end{equation}
where c is just a constant that depends on the bath and will be fixed when we normalize $P(\epsilon)$. Ignorning the higher order terms, we get $P(\epsilon) \propto e^{-\beta\epsilon}$ where we identify $\beta = \frac{\partial\ln\Omega(E)}{\partial E}$.
$\textbf{Question 2}$: Why can we ignore the higher order Taylor terms? I understand that $\epsilon$ is small compared to the heat bath energy $E$ but why am I comparing it to $E$ to get the right measure of the error in probability? 
 A: I think the confusion here has to do with what the Boltzmann distribution describes. It does not give you the probability of finding your small system with a particular energy. Instead, it tells you the probability of finding it in a particular microstate. 
If you want to know the probability of getting a particular energy, you have to sum the Boltzmann probability over the degenerate microstates. This is how you get the Maxwell-Boltzmann distribution.
So, a more correct way to write down what you have above is
$$
P_s \propto \Omega_b\left(E - \epsilon\right)\,,
$$
where $s$ is some microstate of the small system. This leads you to your formula.
Concerning the second order of $\epsilon$. As you say, in equilibrium, we define $\beta = \frac{\partial \ln\Omega\left(E\right)}{\partial E}$. Since $\beta = \frac{1}{k_B T}$ (a constant), the second derivative of the logarithm vanishes, taking away $\epsilon^2$.
A: Question 1
The quantity
    $$P(\epsilon)={\Omega_s(\varepsilon)\Omega_b(E-\varepsilon)\over\Omega_t(\varepsilon)}$$
should be interpreted as the probability of the macroscopic state for which the total energy $E$ is divided in $\varepsilon$ for the system and $E-\varepsilon$ for the bath. Performing a Taylor expansion of the logarithm of the last term, one gets
   $$P(\epsilon)={\rm Cst}\times \Omega_s(\varepsilon)e^{-\beta \epsilon}$$
where $\Omega_s(\varepsilon)$ is the number of microscopic states of the system of energy $\epsilon$. This probability is also the sum of the probabilities of each individual microscopic states whose energies $\epsilon_i$ are equal to $\epsilon$:
   $$P(\epsilon)=\sum_{i/\epsilon_i=\epsilon} p_i$$
Since these microscopic states are equiprobable, $p_i=p(\epsilon)$ for all $i=1,\ldots,\Omega_s(\varepsilon)$ and $P(\epsilon)=\Omega_s(\epsilon)p(\epsilon)$. The probability to find the system in the $i$-th microscopic state is
   $$p_i=p(\epsilon_i)={P(\epsilon)\over\Omega_s(\epsilon)}
   ={\rm Cst}\times e^{-\beta \epsilon_i}$$
I could have used this at the beginning, before performing the Taylor expansion
   $$p_i={\Omega_b(E-\varepsilon_i)\over\Omega_t(\varepsilon_i)}$$
