An approach of successive refinement, based on the Stolz–Cesàro theorem (SCT).
(See N. G. de Bruijn, Asymptotic Methods in Analysis, section $8.6$; also used here.)
Consider (the asymptotics of) $b_n=e^{a_n}$. We have $b_{n+1}=b_n e^{1/b_n}$ or $$b_{n+1}-b_n=b_n(e^{1/b_n}-1).$$
Now $\lim\limits_{n\to\infty}b_n=\infty$ implies $\lim\limits_{n\to\infty}(b_{n+1}-b_n)=1$, hence $\lim\limits_{n\to\infty}b_n/n=1$ by SCT.
So, we put $b_n=n+c_n$ with $c_n=o(n)$. The recurrence becomes $$c_{n+1}-c_n=(n+c_n)(e^{1/(n+c_n)}-1)-1=\frac12\frac1{n+c_n}+O\left(\frac1{n^2}\right),$$ which gives $$\frac12=\lim_{n\to\infty}n(c_{n+1}-c_n)=\lim_{n\to\infty}\frac{c_{n+1}-c_n}{\ln(n+1)-\ln n}=\lim_{n\to\infty}\frac{c_n}{\ln n},$$ the last step again by SCT. Put $c_n=\frac12\ln n+d_n$ with $d_n=o(\ln n)$; then $$d_{n+1}-d_n=\frac1{2n+\ln n+2d_n}-\frac12\ln\frac{n+1}{n}+O\left(\frac1{n^2}\right)=O\left(\frac{\ln n}{n^2}\right),$$ thus $\sum_{n=1}^\infty(d_{n+1}-d_n)$ converges and $\color{red}{\lim\limits_{n\to\infty}d_n=\lambda}$ exists (note that, if we let $a_0$ vary, then $\lambda=\lambda(a_0)$ depends on $a_0$, and this $\lambda$ can be explored in more detail; in particular, we have $\lambda(x+e^{-x})=\lambda(x)+1$ for $x\in\mathbb{R}$).
Continuing, we put $d_n=\lambda+r_n$ and obtain $$r_n-r_{n+1}=\frac14\frac{\ln n}{n^2}+O\left(\frac1{n^2}\right)\implies r_n=\frac14\frac{\ln n}{n}+O\left(\frac1n\right).$$
This process can be developed further, leading to $$\small b_n\asymp n+\frac12\ln n+\lambda+\frac{3\ln n+(6\lambda-2)}{12n}-\frac{3\ln^2n+(12\lambda-10)\ln n+(12\lambda^2-20\lambda+7)}{48n^2}+\dots$$ which, after substituting into $a_n=\ln b_n$, results in
\begin{align}
a_n&\asymp\ln n+\frac{\ln n+2\lambda}{2n}-\frac{3\ln^2n+(12\lambda-6)\ln n+(12\lambda^2-12\lambda+4)}{24n^2}
\\&+{\small\frac{2\ln^3n+(12\lambda-9)\ln^2n+(24\lambda^2-36\lambda+14)\ln n+(16\lambda^3-36\lambda^2+28\lambda-7)}{48n^3}+\dots}
\end{align}
