In Prove $a_t \leq \frac{c}{t}$ I showed that the solutions to $a_0 \in (0, 1]$, $a_n^2 = (1 - a_n)a_{n-1}^2$ satisfy $a_n \le \dfrac{2}{n}$.
Numerical experiments seem to show that $na_n \to 2$ for any initial $a_0 \in (0, 1]$. However, I have not been able to prove this.
These experiments also seem to show that $n(2-na_n)$ does not seem to approach a limit, although it does seem to increase slowly, and this may be due to the limited precision (double precision) used.
So, my questions are
(1) does $na_n \to 2$? (I feel that the answer to this is certainly yes)
(2) Does $\lim_{n \to \infty} n(2-na_n)$ exist? (My feeling is probably)
(3) is there an expansion $a_n =\sum_{k=1}^m \dfrac{c_k}{n^k} $ where $c_1 = 2$?
This may help: My proof that $a_n \le \dfrac{2}{n} $ was by induction based on rewriting the iteration as $a_n =\dfrac{2}{1+\sqrt{1+4/a_{n-1}^2}} $.
I will answer (1) and (2) and try to explain a potential strategy to get the finer results of (3).
First of all let $$f(x) = \frac{2}{1 + \sqrt{1 + \frac{4}{x^2}}},$$ then the sequence $(a_n)$ is defined by $$a_0 > 0, \forall n \geq 0, a_{n+1} = f(a_n).$$
Now you have already shown that $\lim_{n \rightarrow \infty} a_n = 0$, to get finer results we study $f$ at $0$.
Basic computations give :
$$f(x) = x - \frac{1}{2}x^2 - \frac{1}{8}x^3 + o_0(x^3)$$
Then set $b_n = \frac{1}{a_{n+1}} - \frac{1}{a_n}$, using our approximation we have :
\begin{align} \frac{1}{a_{n+1}} &= \frac{1}{a_n - \frac{1}{2}a_n^2 + o(a_n^2)} \\ &= \frac{1}{a_n}(1 + \frac{1}{2}a_n + o(a_n)) \\ &= \frac{1}{a_n} + \frac{1}{2} + o(1). \end{align}
Thus $b_n \rightarrow 2$, and the series $\sum_n b_n$ diverges, using some comparison theorems for divergent series we have :
\begin{align} \sum_{k = 1}^{n} b_k &\sim \frac{n}{2} \\ \frac{1}{a_{n+1}} - \frac{1}{a_0} &\sim \frac{n}{2}. \end{align}
And finally $a_n \sim \frac{2}{n}$ answering (1). Now question (2) is about a second order term, to do this we redo the same computation using the third order approximation of $f$ at $0$. We have :
\begin{align} \frac{1}{a_{n+1}} &= \frac{1}{a_n}\frac{1}{1 - \frac{1}{2}a_n - \frac{1}{8}a_n^2 + o(a_n^2)} \\ &= \frac{1}{a_n} + \frac{1}{2} + \frac{3}{8}a_n + o(a_n). \end{align}
Now the estimate $a_n \sim \frac{2}{n}$ yields $b_n - \frac{1}{2} \sim \frac{3}{8}a_n \sim \frac{3}{4n}$. Using some comparison theorem for divergent series (as before) we get
$$ \frac{1}{a_n} - \frac{n}{2} \sim \frac{3}{4}\ln(n). $$
Finally we need only to inverse this relation, it translates as :
\begin{align} \frac{1}{a_n} - \frac{n}{2} &= \frac{3}{4}\ln(n) + o(\ln(n)) \\ a_n &= \frac{1}{\frac{n}{2} + \frac{3}{4}\ln(n) + o(\ln(n))} \\ &= \frac{2}{n} - 3\frac{\ln(n)}{n} + o(\frac{\ln(n)}{n}). \end{align}
So $n(2 - na_n) \sim 3\ln(n)$ and the limit does not exists. If you want to compute the higher order terms you need to compute the higher order terms of $f$ at $0$ and plug them into the equation, but instead of a theorem for divergent series you will probably need to use the version for convergent series.
The general strategy to determine an approximation of a sequence defined as $u_{n+1} = f(u_n)$, with $f(x) = x - cx^{\alpha} + o_0(x^{\alpha})$ for some $c, \alpha > 0$ is to find a suitable $\gamma$ such that $u_{n+1}^{\gamma} - u_n^{\gamma}$ converges to some nonzero real and then use comparison theorem (or Cesaro's theorem also works in this situation) for convergent/divergent series.