Asymptotic behavior of $a_{n+1}=\frac{a_n^2+1}{2}$ Define a sequence as follows:
$$a_0=0$$
$$a_{n+1}=\frac{a_n^2+1}{2}$$
I would like to know the asymptotic behavior of $a_n$. I already know (by roughly approximating $a_n$ with a differential equation) that
$$a_n\sim 1-\frac{2}{n}$$
as $n\to\infty$. However, my approximation is very crude. Can anyone find a couple more terms? I expect (from numerical data) that the next term is something like $\frac{\log(n)}{n^2}$.
In case anyone wants context for this problem, I am trying to find the optimal strategy for a game with the following rules:

There are $n$ offers of money whose amounts are hidden from you, and whose quantities are random (independently and uniformly distributed from $0$ to $1$). One at a time, the offers are shown to you, and as you view each offer, you may either accept it or reject it. Once you accept an offer, the game is over and you may accept no more offers.

It turns out that $a_n$ is the minimum value of the first offer for which you should accept that offer, in a game with $n+1$ offers. This is why I am interested in the asymptotic behavior of $a_n$.

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}$.

