This answer doesn't contain a closed form expression for sequence $a_n$. I don't believe it has one. In this answer, we try to study how fast $a_n$ decreases to $0$.
Let $f(x) = x e^{-x}$. Notice $0 < f(x) < x$ whenever $x > 0$, the sequence $(a_n)$ defined by the recurrence relation
$$a_n = \begin{cases} 1, & n = 0 \\ f(a_{n-1}), & n > 0\end{cases}$$
is a strictly decreasing sequence bounded from below by $0$. Since $f(x)$is continuous, limit
$a_\infty \stackrel{def}{=}\lim\limits_{n \to \infty} a_n$ exists and $f(a_\infty) = a_\infty \implies a_\infty = 0$. This establish $a_n$ decreases to $0$.
To study how fast $a_n$ decreases to $0$, we look at sequence $\left(a_n^{-1}\right)$ instead. For $n \ge 0$, we have
$$\frac{1}{a_{n+1}} = \frac{1}{a_n} e^{a_{n}} \ge \frac{1}{a_n}(1+a_n) = \frac{1}{a_n} + 1$$
This implies
$$\frac{1}{a_n} = \frac{1}{a_0} + \sum_{k=0}^{n-1} \left( \frac{1}{a_{k+1}} - \frac{1}{a_k}\right) \ge 1 + \sum_{k=0}^{n-1} 1 = 1 + n$$
and hence $a_n$ decreases to $0$ at least as fast as $\frac{1}{n}$.
Let
$$\Delta(x) = \frac{e^x - 1}{x} - 1 - \frac{x}{2} = \sum_{n=2}^\infty \frac{x^n}{(n+1)!}$$
Based on its power series expansion, we know $\Delta(x)$ is strictly increasing, positive and goes to zero as $O(x^2)$ for small $x$. Rewrite the recurrence relation once again, we get
$$\frac{1}{a_{n+1}} = \frac{1}{a_n} e^{a_n} = \frac{1}{a_n} + 1 + \frac{a_n}{2} + \Delta(a_n)
$$
This leads to
$$\begin{align}\frac{1}{a_n} &= n + 1 + \sum_{k=0}^{n-1}\left(\frac{a_k}{2} + \Delta(a_k)\right)
\le n + 1 + \sum_{k=1}^n \left[\frac{1}{2k} + \Delta\left(\frac1k\right)\right]\\
&\le n + 1 + \frac12 H_n + C_0
\end{align}
$$
where $H_n = \sum\limits_{k=1}^n \frac{1}{k}$ is the $n^{th}$ Harmonic number and $$C_0 = \sum_{k=1}^\infty \Delta\left(\frac1k\right) = \sum_{n=2}^\infty \frac{\zeta(n)}{(n+1)!} \approx 0.33493$$
It is known that
$H_n - \log(n+1)$ increases to Euler's constant $\gamma \approx 0.5772156649$ as $n \to \infty$.
This leads to
$$\frac{1}{a_n} \le n + 1 + \frac12 \log(n+1) + C_1
\quad\text{ where }\quad C_1 = \frac{\gamma}{2} + C_0 \approx 0.62354
$$
Plug this upper bound back into the recurrence relation, we obtain
$$\begin{align}\frac{1}{a_{n+1}} - \frac{1}{a_n} - 1 - \frac{1}{2(n+1)}
&= \frac{1}{2a_n} - \frac{1}
{2(n+1)} + \Delta(a_n)\\
&\ge \frac12\left(\frac{1}{n + 1 + \frac12\log(n+1)+C_1} - \frac{1}{n+1}\right)\\
&\ge -\frac14\frac{\log(n+1) + 2C_1}{(n+1)^2}
\end{align}
$$
Summing over different $n$, we can bound $\frac{1}{a_n}$ from both sides
$$-C_2 \le \frac{1}{a_n} - (n + 1 + \frac12 H_n) \le C_0
\quad\text{ where }\quad
C_2 = \frac14\sum_{k=1}^\infty \frac{\log(k) + 2C_1}{k^2} \approx 0.74723$$
If one look at error terms "within" $C_0$ and $C_2$ carefully, all of them have same
signs. The partial sums of these error term converge. As a result, following limit exists
$$C_\infty = \lim_{n\to\infty} \left[ \frac{1}{a_n} - \left( n + 1 + \frac12 H_n\right)\right]
\quad\in \quad(-C_2, C_0)
$$
The sequence $(a_n)$ decreases to $0$ essentially like
$$a_n \approx\frac{1}{n + \frac12 \log n + C + o(1)}
\quad\text{ where }\quad C = 1 + \frac{\gamma}{2} + C_\infty \in (0.54137,1.62354)
$$
As a double check, following are numerical values of $\epsilon_n \stackrel{def}{=} \frac{1}{a_n} - (n + \frac12\log n )$ for selected $n$. It seems to be consistent with above analysis. For the last column, please refer to update below.
$$\begin{array}{r|ll}
n & \epsilon_n & 1/a_n - 1/a_n^{approx} \\
\hline
5 & 1.452256585157472 & 0.000652156218441\\
10 & 1.389752294761182 & 0.000282806941280\\
20 & 1.349378753640906 & 0.000080828265304\\
50 & 1.318816997359306 & 0.000011644667751\\
100 & 1.306359480815289 & 0.000002370371917\\
200 & 1.299203235322892 & 0.000000448836261\\
500 & 1.294298879997029 & 0.000000045923002\\
1000 & 1.292448816475257 & 0.000000007824838\\
2000 & 1.291435425383042 & 0.000000001296712\\
5000 & 1.290768538614429 & 0.000000000128239\\
10000 & 1.290525249116399 & 0.000000000009095\\
\end{array}
$$
Update
To proceed, we need to estimate the value of $C$. It is difficult to extract this value from $a_n$ directly as the convergence of $\epsilon_n$ to $C$ is very slow. As one can see from below, the leading correction of $\epsilon_n$ from $C$ is of the order $\frac{\log n}{n}$.
In order to estimate $C$, we employ the concept of fractional function iterate.
We look for a function $\varphi : [0,\infty) \to \mathbb{R}$ with $\varphi(0) = 1$ and satisfies the functional equation
$$f(\varphi(x)) = \varphi(x+1)$$
It is easy to see $a_n = \varphi(n)$ for all $n \in \mathbb{N}$. As $x$ increases from $0$ to $\infty$, $\varphi(x)$ decreases from $1$ to $0$. If $\psi : [0,1] \to \mathbb{R}$ is the inverse function of $\varphi$, then $\psi(1) = 0$ and satisfy the Abel equation:
$$\psi(f(z)) = \psi(z) + 1$$
Since $\displaystyle\;\frac{1}{a_n} \sim n + \frac12\log n + C$ for large $n$, we expect
$\psi(z)$ has the form
$$\psi(z) = \frac1z + \frac12 \log z - C + F(z)$$
for some function $F(z)$ which tends to $0$ as $z \to \infty$. Plug this into Abel equation and simplifies, we get
$$\frac1z e^z + \frac12(\log(z) - z) + F(ze^{-z}) = \frac1z + \frac12 \log z + F(z) + 1\\
\iff F(z) - F(ze^{-z}) = \frac1z\left( e^z - 1 - z - \frac{z^2}{2}\right)$$
Notice the cancellation of the problematic $\log z$ factor. There is nothing to stop $F(z)$ from having a power series expansion at $z = 0$. With help of a CAS,
we find
$$\begin{align} F(z) =
\frac{z}{6}
&+\frac{z^2}{16}+\frac{19 z^3}{540}+\frac{z^4}{48}+\frac{41 z^5}{4200}+\frac{37z^6}{103680}-\frac{18349 z^7}{3175200}\\
&-\frac{443 z^8}{80640}+\frac{55721 z^9}{21555072}+\frac{84317 z^{10}}{6912000}
+\frac{2594833561\,{z}^{11}}{399567168000}\\
&-\frac{152043613 z^{12}}{5748019200}-\frac{830066563 z^{13}}{17346700800}+\frac{209959420499 z^{14}}{3905587445760}\\& \;\;\vdots
\end{align}
$$
In principle, we can extract $C$ using the condition $\psi(1) = 0 \iff C = F(1)$.
However, this is not that accurate because the coefficients of $z^k$ in $F(z)$ doesn't decade quickly. A more accurate way is use the relation $\psi(a_n) = n$. This leads to
$$C = \frac{1}{a_n} + \frac12\log(a_n) + F(a_n) - n$$
With $n$ as small as $10$ and keeping $F(z)$ up to $z^{10}$, we already obtain a decent estimate $$C \approx 1.290247208687763$$
Once we have a reasonable estimate of $C$, we can invert the equation
$$n = \psi(z) = \frac{1}{z} + \frac12\log(z) - C + F(z)$$
and substitute the resulting $z$ by $a_n$. Using a CAS again, we find
$$\begin{align}
\frac{1}{a_n} = n &+ \frac12\log n + C +
\frac{1}{12n}\left(6C-2 + 3\log n\right)\\
& - \frac{1}{48n^2}\left(
(12C^2 - 20C + 7) + (12C-10)\log(n) + 3\log(n)^2\right)\\
& \;\;\vdots
\end{align}
$$
Let $a_n^{approx}$ be this approximation of $a_n$, the last column in the table above is the difference between $\frac{1}{a_n}$ and $\frac{1}{a_n^{approx}}$. As one can see, the differences decrease quickly. This partially validate the correctness of this approach and accuracy in estimation of $C$.
Best Answer
This is an answer to part C.
Currently, we know how to iterate two types of quadratic sequences. They are for quadratics in one of these two forms: $$q_1(x)=ax^2+bx+\frac{b^2-2b}{4a}$$ $$q_2(x)=ax^2+bx+\frac{b^2-2b-8}{4a}$$ See the wikipedia page on iterated functions or this blog post for a more in-depth explanation of how to iterate these two types.