Alguien en las redes sociales (donde enlacé este post), Prof. Antoine Chambert-Loir dio una solución. El esquema que sigue procede de él; los errores, si los hay, son míos.
La idea es establecer $$a_n = \operatorname{sl}^2(u_n) \tag{1} $$ donde $\operatorname{sl}$ es el función seno lemniscata .
Entonces la relación de recurrencia pasa a ser $$ \operatorname{sl}^2(u_{n+1}) = \frac{1-\sqrt{1-\operatorname{sl}^2(u_n)}}{1+\sqrt{1+\operatorname{sl}^2(u_n)}} = \operatorname{sl}^2(\tfrac{1}{2}u_{n}) \tag{2} $$ utilizando el fórmulas de bisección e identidades relativas $\operatorname{sl}$ y $\operatorname{cl}$ . Esto demuestra que $u_{n+1} = \frac{1}{2}u_n$ y, por lo tanto $$ u_n = \frac{2u_1}{2^{n}} \tag{3} $$ Pero $a_1=1$ implica $\operatorname{sl}(u_n) =1$ que entonces da $$ u_1 = \int_0^1 \frac{dx}{\sqrt{1-x^4}} = \frac{\sqrt{\pi}\Gamma(5/4)}{\Gamma(3/4)} \tag{4} $$ que finalmente conduce a $$ a_n = \operatorname{sl}^2(u_n) \operatorname*{\sim}_{n\to\infty} u_n^2 = \frac{4u_1^2}{4^n} = \frac{4\pi\Gamma(5/4)^2}{\Gamma(3/4)^2} \cdot 4^{-n} \tag{5} $$ donde $\varpi^2=\frac{4\pi\Gamma(5/4)^2}{\Gamma(3/4)^2}\approx 6.875185818$ ( $\varpi$ siendo el constante lemniscata ).
Lo anterior, según la pregunta del OP, corresponde a la condición inicial $a_1=1$ . Utilizando el mismo argumento, se obtiene fácilmente que, para $a_1=\alpha\in[0,1]$ la constante límite para $4^n a_n$ viene dada por $$ C(\alpha) = 4\operatorname{arcsl}(\sqrt{\alpha})^2 = 4\left(\int_0^{\sqrt{\alpha}} \frac{dx}{\sqrt{1-x^4}}\right)^2 \tag{6} $$ donde $\operatorname{arcsl}$ es el función arcoseno lemniscata . El comportamiento de $C(\alpha)$ como $\alpha$ oscila entre $0$ a $1$ se representa a continuación.