Simple proof that the Airy function decays faster than any polynomial

I will assume that what we know about $y = \operatorname{Ai}(x)$ is that $y'' = xy$ and $ y \to 0$ through positive values as $x \to +\infty$.

Assume $y > 0$ for $x > M$. Let $k$ be fixed and large enough so that $k^2 > M$. Let $y_1 = Ce^{-kx}$, where $C$ is chosen so that $y_1(k^2) = y(k^2)$. We will prove that $y \leq y_1$ for $x \geq k^2$. This shows that $y$ is negligible compared to any exponential function.

The function $y - y_1$ takes the value $0$ at $k^2$ and tends to $0$ as $x \to +\infty$.

Assume for a contradiction that $y > y_1$ for some $x > k^2$. Then $y - y_1$ attains a positive maximum, say at $x_0 > k^2$. Therefore the second derivative of $y - y_1$ must be $\leq 0$ at $x_0$. Hence we obtain the absurd consequence that $$0 \geq y''(x_0) - y_1''(x_0) = x_0 y(x_0) - k^2 y_1(x_0) \geq k^2 [y(x_0) - y_1(x_0)] > 0.$$

This is somewhat weaker than what you want, but may be useful.

To study what happens as $x \to +\infty$, it's convenient to use the change of variables $t = 1/x$ so that $t \to 0+$. The Airy equation becomes

$$ (t^2 u')' = u/t^3 $$

Now suppose a solution $u$ has $u \sim c t^m$ as $t \to 0+$ for some $m, c > 0$. Then $(t^2 u')' \sim c t^{m-3}$.

Note that $t^2 u'$ is increasing. If $u'(t_0) = 0$ for some $t_0 > 0$, we'd have $u'(t_1) < 0$ for $t_1$ with $0 < t_1 < t_0$, and then there would be $\epsilon > 0$ with $u' < - \epsilon/t^2 $ for $0 < t < t_1$. But since $1/t^2$ is not integrable on that interval, this is impossible. Similarly we can't have all $t^2 u' > \epsilon > 0$. So we must have $u' > 0$, and $t^2 u' \to 0$ as $t \to 0+$, and $$ t^2 u'(t) = \int_0^t (s^2 u'(s))' ds \sim c \int_0^t s^{m-3}\; ds$$ In particular $s^{m-3}$ must be integrable on $(0,t]$ so $m > 2$, and $$t^2 u'(t) \sim \dfrac{c}{m-2} t^{m-2} $$
That makes $u'(t) \sim \dfrac{c}{m-2} t^{m-4} $. Again, $t^{m-4}$ must be integrable, so $m > 3$, and $$ u(t) \sim \dfrac{c}{(m-2)(m-3)} t^{m-3} $$ But that contradicts the assumption $u \sim c t^m$.

I haven't thought through a fully rigorous argument, but I think we can get the estimate $$\textrm{Ai}(x) \le \exp\left( - \tfrac23 x^{3/2} +C\sqrt{x}\right) \qquad (x>0)$$ (which holds for all positive $x$, with $C$ being some positive constant) using some fairly low-tech and intuitive bounds. This comes pretty close to the more precise asymptotic formula $$ \textrm{Ai}(x) \sim \frac{1}{2\sqrt{\pi} x^{1/4}} \exp\left( - \tfrac23 x^{3/2} \right) \qquad (x\to\infty) $$ and explains where the exponential decay in $x^{3/2}$ comes from (which is a bit more than what you asked for -- it's certainly faster than any negative power of $x$, and even a bit faster than a simple exponential decay in $x$ because of the $3/2$ exponent), and even gives the correct coefficient $2/3$ preceding $x^{3/2}$.

The idea is to compare $\textrm{Ai}(x)$ to the function $G(x)$ that solves a discretized version of the Airy differential equation, where we replace $$ y''(x) = x y(x) $$ by $$ y''(x) = (n-1) y(x) \qquad \textrm{on any interval }n-1<x<n $$ for integer $n\ge 1$. Specifically I construct this $G(x)$ by gluing together a bunch of exponentials of the form $G(x)=c_n e^{-\sqrt{n-1}x}$ on each interval $[n-1,n]$. The constants $c_n$ are computed inductively from the requirement that $G(x)$ should be continuous, and that we have $G(0)=\textrm{Ai}(0)$. This leads to the piecewise formula $$ G(x) = \begin{cases} \textrm{Ai}(0) & \textrm{if }0\le x<1, \\ \textrm{Ai}(0)\exp\big(-\sqrt{1}(x-1)\big) & \textrm{if }1\le x<2, \\ \textrm{Ai}(0)\exp\big(-\sqrt{1}-\sqrt{2}(x-2)\big) & \textrm{if }2\le x<3, \\ \textrm{Ai}(0)\exp\big(-\sqrt{1}-\sqrt{2}-\sqrt{3}(x-3)\big) & \textrm{if }3\le x<4, \\ \textrm{Ai}(0)\exp\big(-\sqrt{1}-\sqrt{2}-\sqrt{3}-\sqrt{4}(x-4)\big) & \textrm{if }4\le x<5, \\ \cdots \\ \textrm{Ai}(0) \exp\left(-\sum_{k=1}^{n-1} \sqrt{k} - \sqrt{n}(x-n)\right) & \textrm{if }n\le x<n+1, \\ \cdots \end{cases} $$

Here is a plot showing the two functions $\textrm{Ai}(x)$ and $G(x)$:

Now the trick is to observe that $0<\textrm{Ai}(x)\le G(x)$ for any $x>0$. This shouldn't be too hard to prove rigorously; the intuition behind it is that the change from the original Airy equation to the discretized version makes it harder for the function to decay by speeding up the second derivative of $y(x)$ just a bit. (Edit: on second thought, this version of the equation slows down the second derivative, so the comparison is a bit more subtle than I thought at first. Nonetheless, the plot suggests the inequality is true, and if it is then I suspect it should not be hard to justify rigorously.)

Assuming this, we see that, for $x$ in the interval $[n,n+1]$ with $n\ge1$ integer,

$$ \textrm{Ai}(x) \le \textrm{Ai}(0)\exp\left( - \sum_{k=1}^{n-1} \sqrt{k} \right),$$ and this is easily seen to behave like $\exp\left( - \tfrac32 x^{3/2} +C\sqrt{x}\right)$, e.g., by a comparison of the sum inside the exponential with the integral $\int_0^n \sqrt{x}\,dx = \tfrac23 n^{3/2}$.