After doing integration by parts, he ended up with
$$
\operatorname{Ei}(z)=\frac{e^z}{z} \biggl [\sum_{k=0}^n \frac{k!}{z^k}+\epsilon _n(z) \biggr ].
$$
[...] He then showed that $$|\epsilon _n(z)| \leq \frac{(n+1)!}{|z|^{n+1} (\sin \delta)^{n+2}},$$
where $\delta$ is some constant such that it ranges $(0,\pi)$, and such that $|\arg(-z)| \leq \delta - \pi$.
[...] I think I should point out that the $\sin \delta$ above is not in absolute values, which I think is problematic, since provided the range of $\delta $, $\sin \delta $ could take on negative values and the inequality above will not be true anymore.
If $\delta \in (0,\pi)$ then $\sin \delta > 0$, so there's no need to worry about the absolute value.
Does the last result above show that $\epsilon _n (z) = O(z^{-(n+1)})$?
The number $\delta$ is constant with respect to $z$ here, so there really is a constant
$$
C = \frac{(n+1)!}{(\sin \delta)^{n+2}}
$$
such that
$$
|e_n(z)| \leq C |z|^{-(n+1)}.
$$
This is exactly what it means to write $e_n(z) = O(z^{-(n+1)})$.
Now, owing to the fact that there's still that $\frac{e^z}{z}$ factor outside the RHS of the second equation from above, then from what I understood, it should say that:
$$
Ei(z) = \frac{e^z}{z} \sum_{k=0}^n \frac{k!}{z^k} + O \left (\frac{e^z}{z^{(n+2)}} \right),
$$
which doesn't really fit the Poincare asymptotic expansion. Or is because I shouldn't see it in light of the Poincare asymptotic expansion but in terms of some other asymptotic scale?
There are two ways to think about it. One way is to agree that it shows that
$$
\sum_{k=0}^n \frac{k!}{z^k}+\epsilon _n(z) = \sum_{k=0}^n \frac{k!}{z^k}+O(z^{-(n+1)})
$$
is a Poincaré expansion, so what the formula
$$
\operatorname{Ei}(x) \sim \frac{e^x}{x} \sum_{k=0}^{\infty} \frac{k!}{z^k}
$$
means is precisely that
$$
\operatorname{Ei}(x) = \frac{e^x}{x} \left[\sum_{k=0}^{n} \frac{k!}{z^k}+O(z^{-(n+1)})\right]
$$
for every $n \in \mathbb N$.
The other way is to expand things as you've done and write
$$
\operatorname{Ei}(z) = \sum_{k=0}^n \frac{k!e^z}{z^{k+1}} + O\left( \frac{e^z}{z^{(n+2)}} \right). \tag{$*$}
$$
This is now a generalized asymptotic expansion $\sum_k g_k(z)$, where the summands
$$
g_k(z) = \frac{k!e^z}{z^{k+1}}
$$
form an asymptotic sequence. Indeed we have
$$
O\left( \frac{e^z}{z^{(n+2)}} \right) \subset o(g_n(z))
$$
as $z \to \infty$, so a consequence of $(*)$ is that
$$
\operatorname{Ei}(z) = \sum_{k=0}^n g_k(z) + o(g_n(z))
$$
for every fixed $n \in \mathbb N$ as $z \to \infty$. This is precisely the statement that
$$
\operatorname{Ei}(z) \sim \sum_{k=0}^\infty g_k(z)
$$
as $x \to \infty$, that is
$$
\operatorname{Ei}(z) \sim \sum_{k=0}^\infty \frac{k!e^z}{z^{k+1}} \tag{$**$}
$$
as $x \to \infty$. Now, the interpretation would be equivalent either way (i.e. the Big O statements one could obtain from the asymptotic series) so one could formally factor out $e^z/z$ from this "sum" and thus write the result as
$$
\operatorname{Ei}(z) \sim \frac{e^z}{z}\sum_{k=0}^\infty \frac{k!}{z^k} \tag{$***$}
$$
as $z \to \infty$.
You should convince yourself that $(**)$ and $(***)$ really are equivalent statements.
Best Answer
As $k>n-1$ the integral in the decomposition you proposes diverges. You can write the original integral as $$\exp(-x)\int_0^\infty (1+u)^{-n}\exp(-ux)du$$ and integrate by parts several times. Firstly $$E_n(x)=e^{-x}\left[\frac{1}{x}-\frac{n}{x}\int_0^\infty (1+u)^{-n-1}e^{-ux}du\right]$$ which gives the first order term. Going on gives you the following terms. Up to the first order $$E_n(x)=\frac{e^{-x}}{x}+O(x^{-2}e^{-x})$$