Your asymptotic approximation is correct :
$$y (x)\sim C_2 (ax)^{-1/4} \exp \left[-\frac{1}{\epsilon} \sqrt{a} \frac{2}{3} x^{3/2}\right] \tag 1$$
according to the boundary condition $y(\infty)=0$.
Note : A second order ODE requires at least two boundary conditions to expect having a unique solution. Since the wording of the question specifies only the condition $y(\infty)=0$, the solution remains undetermined : The coefficient $C_2$ cannot be computed.
Another key point is that an asymptotic formula is valid only for $x\to\infty$, not for $x\to 0$.So, don't expect Eq.$(1)$ makes sens at $x=0$. It is normal that the asymptotic approximation be far to be correct for small values of $x$.
If you want an approximate around $x=0$, expand $y(x)$ in power series.
$$y(x)\simeq c_0+c_1x+c_2x^2+...$$
Of course, this formula is not valid for large $x$, so the boundary condition at infinity cannot be fulfilled. Moreover, two boundary conditions should be required close to $x=0$ in order to determine the coefficients $c_0$ , $c_1$, $c_2$ , …
In fact, as Claude Leibovici quite rightly wrote in comments, the exact solving involves the Airy functions.
$$y(x)=c_1\text{Ai}\left(\sqrt[3]{\frac{a}{\epsilon^2}}\:x\right)+c_2\text{Bi}\left(\sqrt[3]{\frac{a}{\epsilon^2}}\:x\right)$$
Ai and Bi are the Airy functions http://mathworld.wolfram.com/AiryFunctions.html
The condition $y(\infty)=0$ is satisfied any $c_1$ and $c_2$, which shows that the condition $y(\infty)=0$ is not sufficient to fully determine a solution.
On physicists viewpoint, defining the properties of a boundary layer is certainly a way to introduce some boundary conditions in the range around $x=0$. Supposing that $y(x)\simeq y_0+y'_0 x$ in a layer of a given thickness $\delta$ allows to compute $c_1$ and $c_2$ as a function of $y_0$ , $y'_0$ and $\delta$, thanks to the series expansion of the Airy functions. Then the asymptotic expansion of the Airy functions will leads to the coefficient $C_2$ in the asymptotic approximation $(1)$.
http://functions.wolfram.com/Bessel-TypeFunctions/AiryAi/06/
http://functions.wolfram.com/Bessel-TypeFunctions/AiryBi/06/
For the first question, the answer is yes. For the third question (i.e. the bold question), yes it works because you get $\frac{d(n)}{n^{\varepsilon}}\leqslant C_{\varepsilon}$ where $\displaystyle C_{\varepsilon}=\prod_{p<e^{1/\varepsilon}}C_{p,\varepsilon} $ does not depend of $n$. As for $2.$, I'll use $3.$. First notice that the bound $C_{\varepsilon}:=\sup\limits_{n\geqslant 1}\frac{d(n)}{n^{\varepsilon}}$ is reached for $\displaystyle n_{\varepsilon}:=\prod_p p^{\alpha_{p,\varepsilon}}$ where $\displaystyle\alpha_{p,\varepsilon}:=\left\lfloor \frac{1}{p^{\varepsilon}-1} \right\rfloor$. Indeed, if $f(\alpha):=\frac{\alpha+1}{p^{\alpha\varepsilon}}$, then
$$ f(\alpha+1)\geqslant f(\alpha)\iff\frac{\alpha+2}{\alpha+1}\geqslant p^{\varepsilon}\iff\alpha\leqslant\alpha_{p,\varepsilon} $$
and $f$ reaches its maximum at $\alpha=\alpha_{p,\varepsilon}$. Now let $x_k:=\left(1+\frac{1}{k}\right)^{1/\varepsilon}$, then
$$ \alpha_{p,\varepsilon}=k\iff \frac{1}{p^{\varepsilon}-1}-1<k\leqslant\frac{1}{p^{\varepsilon}-1}\iff x_{k+1}<p\leqslant x_k $$
for $k\geqslant 1$. Let $k_0:=\alpha_{2,\varepsilon}$, then there is no prime $p$ such that $p\leqslant x_{k_0+1}$ because $x_{k_0+1}<2$ so we have
$$ n_{\varepsilon}=\prod_{k\leqslant k_0}\left(\prod_{x_{k+1}<p\leqslant x_k}p\right)^k $$
From this expression we can deduce the two following estimations :
$$ \ln n_{\varepsilon}=\vartheta(x_1)+\mathcal{O}\left(x_1^{3/4}\right) \ \ \text{ and }\ \ \ln d(n_{\varepsilon})=(\ln 2)\pi(x_1)+\mathcal{O}\left(x_1^{3/4}\right) $$
Indeed, $\displaystyle\ln n_{\varepsilon}=\sum_{k\leqslant k_0}k(\vartheta(x_k)-\vartheta(x_{k+1}))=\sum_{k\leqslant k_0}\vartheta(x_k)$ and, using $x_2\leqslant x_1^{\frac{\ln 3}{\ln 2}-1}$, we have
$$ \sum_{2\leqslant k\leqslant k_0}\vartheta(x_k)\ll k_0\vartheta(x_2)\ll k_0 x_2\ln x_2\ll x_2(\ln x_2)^2\ll x_1^{\frac{\ln 3}{\ln 2}-1}(\ln x_1)^2\ll x_1^{3/4} $$
because $\frac{\ln 3}{\ln 2}-1\approx 0.58\leqslant 0.75$. As for the other approximation, we have
$$ \ln d(n_{\varepsilon})=\sum_{k\leqslant k_0}\ln(k+1)(\pi(x_k)-\pi(x_{k+1}))=\sum_{k\leqslant k_0}\ln\left(1+\frac{1}{k}\right)\pi(x_k) $$
and
$$ \sum_{2\leqslant k\leqslant k_0}\ln\left(1+\frac{1}{k}\right)\pi(x_k)\ll\frac{k_0 x_2}{\ln x_2}\ll x_1^{3/4} $$
using the same arguments. Now let $R(x)$ be such that $\pi(x)-{\rm li}(x)\ll R(x)$ and $\vartheta(x)-x\ll R(x)$, then
$$ \ln d(n)\leqslant \ln C_{\varepsilon}+\varepsilon\ln n\leqslant\ln d(n_{\varepsilon})-\varepsilon\ln n_{\varepsilon}+\varepsilon\ln n\leqslant(\ln 2)\pi(x_1)-\varepsilon\vartheta(x_1)+\varepsilon\ln n+\mathcal{O}\left(x_1^{3/4}\right)$$
We then use the bound $R(x)\gg x^{4/5}$ and we obtain
$$ \ln d(n)\leqslant (\ln 2){\rm li}(x_1)-\varepsilon x_1+\varepsilon\ln n+\mathcal{O}(R(x)) $$
In order to cancel the first order terms, we chose $\varepsilon:=\frac{\ln 2}{\ln\ln n}$ so that $x_1=2^{1/\varepsilon}=\ln n$. We thus have $\ln d(n)\leqslant(\ln 2){\rm li}(\ln n)+\mathcal{O}(R(\ln n))$ and, using the bound $R(x)\ll\frac{x}{(\ln x)^2}$, we finally get
$$ \ln d(n)\leqslant \frac{(\ln 2)\ln n}{\ln \ln n}+\mathcal{O}\left(\frac{\ln n}{(\ln \ln n)^2}\right)\underset{n\rightarrow +\infty}{\sim}\frac{(\ln 2)\ln n}{\ln \ln n} $$
Best Answer
Let's look at them in turn.
We want to show that $$ O\left(\frac{1}{\epsilon}\right)^{\exp(1/\epsilon)}=\exp(\exp(O(1/\epsilon))). $$
For ease of notation, I drop most $O(\cdot)$ constants. Then note that $$ \left( \frac{1}{e} \right) ^{ \exp(1/\epsilon) } = e^{\log(1/\epsilon)\cdot \exp(1/\epsilon)}$$
and in general, $\log x \cdot e^x = O(e^{x +\delta})$ for any $\delta > 0$. Thus $$ \left( \frac{1}{e} \right) ^{ \exp(1/\epsilon) } = e^{\log(1/\epsilon)\cdot \exp(1/\epsilon)} = \exp (\, \exp ( \, O (1/\epsilon )\, )\, )$$
Let's follow his suggestion. So set $\epsilon = C/\log\log n$ for some large $C$. Then we want to compare the growth of $\exp ( \exp ( O(1/\epsilon)))$ and $n^\epsilon$. I again drop most $O(\cdot)$ constants for ease. The first one:
$$ \exp( \exp ( 1/\epsilon)) = \exp ( \exp ( \log \log n / C) ) = \exp \exp (\log (\log n)^{1/C} ) = \exp ( (\log n)^{1/C}).$$
The second one stays:
$$ n^{C/\log \log n}$$
Now take logs. Then we compare
$$ (\log n)^{1/C} \quad \text{and} \quad \frac{C}{\log \log n} \log n$$
Similar to the fact above, $\frac{x}{\log x} > x^{1 - \delta}$ for any $\delta > 0$. So for large $C$, the left side is like a small power (much less than $1$, say) of $\log n$ and the right side is like a larger power (pretty close to $1$) of $\log n$.
So that's one way of getting his two results. $\diamondsuit$