In the present case, there is a shortcut to Watson's lemma (mentioned in the comments), which is to scale properly the variable of integration of the integral $I(x)$ to be evaluated.
Let $x=1/z^2$ with $z\gt0$, hence $z\to0^+$. Using the change of variable $t\to z^2t$, one gets
$$
I(x)=z^2\int_0^{+\infty}\mathrm e^{-t}\log(1+z\sqrt{t})\,\mathrm dt.
$$
For every $N\geqslant0$, an expansion of $\log(1+s)$ up to order $s^N$ when $s\to0$ is
$$
\log(1+s)=\sum_{n=1}^N(-1)^{n-1}\frac1ns^n+o(s^{N}).
$$
This yields
$$
I(x)=z^2\sum_{n=1}^N(-1)^{n-1}I_nz^n+o(z^{N+2})=\sum_{n=1}^N(-1)^{n-1}I_nx^{-1-n/2}+o(x^{-1-N/2}),
$$
where, for every $n\geqslant1$,
$$
I_n=\frac1n\int_0^{+\infty}\mathrm e^{-t}t^{n/2}\,\mathrm dt=\frac1n\Gamma\left(\frac{n}2+1\right)=\frac12\Gamma\left(\frac{n}2\right).
$$
Finally, for every $N\geqslant0$,
$$
I(x)=\sum_{n=1}^N(-1)^{n-1}\frac12\Gamma\left(\frac{n}2\right)x^{-1-n/2}+o(x^{-1-N/2}).
$$
Note that the radius of convergence of the series $\sum\limits_n\Gamma\left(\frac{n}2\right)t^n$ being zero, the formulas above are indeed asymptotic expansions.
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
Here I will solve the integral in terms of complete elliptic integrals. Then the claim follows from the known asymptotic behaviour of these special functions.
Note that the integrand and the bounds are symmetric about $x = 1/2$. This becomes even clearer if we apply the shift $x = y + 1/2$. Writing $\delta = \varepsilon/2$ temporarily to simplify notation a bit, we find \begin{align} I (\varepsilon/2) &= \int \limits_{-(1-\varepsilon)/2}^{(1-\varepsilon)/2} \mathrm{d} y \, \frac{1 - 4 y^2}{\sqrt{[(1+\varepsilon)^2 - 4 y^2][(1-\varepsilon)^2 - 4 y^2]}} \\ &= 2 \int \limits_0^{(1-\varepsilon)/2} \mathrm{d} y \, \frac{1 - 4 y^2}{\sqrt{[(1+\varepsilon)^2 - 4 y^2][(1-\varepsilon)^2 - 4 y^2]}}\, . \end{align} Rescaling with $y = (1-\varepsilon)z/2$ and rearranging, we end up with \begin{align} I (\varepsilon/2) &= \int \limits_0^1 \mathrm{d} z \, \frac{1 - (1-\varepsilon)^2 z^2}{\sqrt{[(1+\varepsilon)^2 - (1-\varepsilon)^2 z^2][1-z^2]}} \\ &= \int \limits_0^1 \mathrm{d} z \, \left(\frac{(1+\varepsilon) \sqrt{1 - \left(\frac{1-\varepsilon}{1+\varepsilon}\right)^2 z^2}}{\sqrt{1-z^2}} - \frac{\varepsilon (2 + \varepsilon)}{(1+\varepsilon)\sqrt{1 - \left(\frac{1-\varepsilon}{1+\varepsilon}\right)^2 z^2}\sqrt{1-z^2}} \right)\, . \end{align} But these are just two complete elliptic integrals, which we can simplify using Landen transformations: \begin{align} I (\varepsilon/2) &= (1+\varepsilon) \operatorname{E} \left(\frac{1-\varepsilon}{1+\varepsilon}\right) - \frac{\varepsilon (2 + \varepsilon)}{1 + \varepsilon} \operatorname{K} \left(\frac{1-\varepsilon}{1+\varepsilon}\right) \\ &= \operatorname{E} \left(\sqrt{1-\varepsilon^2}\right) - \frac{\varepsilon^2}{2} \operatorname{K} \left(\sqrt{1-\varepsilon^2}\right) \, . \end{align}
Finally, their asymptotic expansions yield \begin{align} I(\delta) &= \operatorname{E} \left(\sqrt{1-4\delta^2}\right) - 2 \delta^2 \operatorname{K} \left(\sqrt{1-4\delta^2}\right) \\ &\sim 1 - \delta^2 - \left(\log(\delta/2) + \frac{5}{4}\right) \delta^4 + \operatorname{\mathcal{O}} \left(\log(\delta) \delta^6\right) \end{align} as $\delta \to 0^+$.