Write $a_n = 1 - \epsilon_n$ and notice that $(\epsilon_n)$ solves
$$ \epsilon_{n+1} = \epsilon_n - \frac{\epsilon_n^2}{2}.$$
For the purpose of future use, we allow $\epsilon_0$ to take any value in $(0, 1]$. This type of sequence is well-studied, and here is a method of extracting asymptotic forms up to certain order in a bootstrapping manner.
Since $0 \leq \epsilon_n \leq 1$ and $(\epsilon_n)$ is monotone decreasing, we have $\epsilon_n \to 0$ as $n\to\infty$.
We have $ \frac{1}{\epsilon_{n+1}} - \frac{1}{\epsilon_n} = \frac{1}{2-\epsilon_n} $. So by Stolz–Cesàro theorem (a.k.a. L'Hospital's theorem for sequence),
$$ \lim_{n\to\infty} \frac{1/\epsilon_n}{n} = \lim_{n\to\infty} \left(\frac{1}{\epsilon_{n+1}} - \frac{1}{\epsilon_n} \right) = 2.$$
This shows that
$$\epsilon_n = (1 + o(1))\frac{2}{n}. $$
Pushing this idea further, we may utilize $ \frac{1}{\epsilon_{n+1}} - \frac{1}{\epsilon_n} = \frac{1}{2} + \frac{\epsilon_n}{2(2-\epsilon_n)} $. Again, by Stolz–Cesàro theorem,
$$ \lim_{n\to\infty} \frac{\frac{1}{\epsilon_n} - \frac{n}{2}}{\log n}
= \lim_{n\to\infty} \frac{\frac{1}{\epsilon_{n+1}} - \frac{1}{\epsilon_n} - \frac{1}{2}}{\log (n+1) - \log n}
= \frac{1}{2}, $$
and so, $\frac{1}{\epsilon_n} = \frac{n}{2} + \frac{1 + o(1)}{2}\log n$. This in turn implies
$$ \epsilon_n = \frac{2}{n} - (2 + o(1)) \frac{\log n}{n^2}. $$
Finally, by writing
\begin{align*}
\frac{1}{\epsilon_{n+1}} - \frac{1}{\epsilon_n}
&= \frac{1}{2} + \frac{1}{2(n+1)} + \underbrace{\left( \frac{\epsilon_n}{2(2-\epsilon_n)} - \frac{1}{2(n+1)} \right)}_{=\mathcal{O}(n^{-2})}, \\
\end{align*}
we obtain
\begin{align*}\frac{1}{\epsilon_n}
&= \frac{1}{\epsilon_0} + \frac{n}{2} + \frac{H_n}{2} + \sum_{k=0}^{n-1} \left( \frac{\epsilon_k}{2(2-\epsilon_k)} - \frac{1}{2(k+1)} \right) \\
&= \frac{n}{2} + \frac{H_n}{2} + C\left(\epsilon_0\right) + \mathcal{O}_{\epsilon_0}\left(\frac{1}{n}\right). \tag{*}
\end{align*}
Here, $H_n = \sum_{k=1}^{n} \frac{1}{k}$ is the harmonic number and $C(\epsilon_0)$ is a constant depending only on $\epsilon_0$. Also, the implicit bound of $\mathcal{O}_{\epsilon_0}$ depends only on $\epsilon_0$. Moreover, $C(\cdot)$ admits a nice property. Let $f(x) = x-x^2/2$ so that $\epsilon_{n+1} = f(\epsilon_n)$. Then
$$ C(x) = \frac{1}{\epsilon_0} + \sum_{k=0}^{\infty} \left( \frac{f^{\circ k}(x)}{2(2-f^{\circ k}(x))} - \frac{1}{2(k+1)} \right), $$
where $f^{\circ k}$ denotes the $k$-fold function composition of $f$ together with $f^{\circ 0} = \mathrm{id}$. Then $\text{(*)}$ can be recast to
$$ \frac{1}{f^{\circ n}(x)} = \frac{n}{2} + \frac{H_n}{2} + C(x) + \mathcal{O}_{x}\left(\frac{1}{n}\right). $$
Then it follows that
\begin{align*}
0
&= \frac{1}{f^{\circ (n+1)}(x)} - \frac{1}{f^{\circ (n+1)}(x)} \\
&= \left( \frac{n+1}{2} + \frac{H_{n+1}}{2} + C(x) + \mathcal{O}\left(\frac{1}{n}\right) \right) - \left( \frac{n}{2} + \frac{H_n}{2} + C(f(x)) + \mathcal{O}\left(\frac{1}{n}\right) \right),
\end{align*}
and letting $n\to\infty$ gives
$$ C(f(x)) = C(x) + \frac{1}{2}. $$
Summarizing,
Let $(\epsilon_n)$ be a sequence which solves $\epsilon_{n+1} = f(\epsilon_n)$, where $\epsilon_0 \in (0, 1]$ and $f(x) = x - x^2/2$. Then there exists a function $C : (0, 1] \to \mathbb{R}$ such that
$$ \frac{1}{\epsilon_n} = \frac{n}{2} + \frac{H_n}{2} + C(x) + \mathcal{O}_x\left(\frac{1}{n}\right), $$
and so,
$$ \epsilon_n = \frac{2}{n} - \frac{2H_n}{n^2} - \frac{4C(x)}{n^2} + \mathcal{O}_x\left( \frac{\log^2 n}{n^3} \right). $$
Moreover, $C$ solves $C(f(x)) = C(x) + \frac{1}{2}$.