With regard to the first question:
It's not so much the rapid oscillations themselves, but rather that the area associated with each oscillation decreases as $\lvert x\rvert \rightarrow \infty $.
Considering the real part of the integral:
$$ \Re(I) = \int_{-\infty}^{\infty} \cos{(x^3/3+sx)}dx = 2\int_{0}^{\infty} \cos{(x^3/3+sx)}dx$$
This integral can be realized as a summation of a sequence $(a_n)$, where $a_n$ is the area of each half-oscillation:
$$ \Re(I)=2\sum_{n=0}^{\infty} a_n$$
For this summation:
- $a_n$ will alternate between positive and negative values
- You can show that $a_n \rightarrow 0$
- $\lvert a_n \rvert$ is decreasing
And hence by the alternating series test, the summation converges
This material can be found in Bender and Orszag, Chapter 3. (They use the probabilist's convention, though, and derive $He(x)$ rather than $H(x)$.)
There is an alternative way to go about analyzing this problem in which we compute an asymptotic series about the irregular point at infinity. This first involves extracting the leading order behavior of the solutions to the original (scaled) differential equation, given by
$$
-\frac{1}{2}y''+\frac{1}{2}x^2y-\epsilon y=0\,,
$$
by first making the transformation $y(x)=e^{S(x)}$ and then using the method of dominant balance. This method is systematic if not entirely rigorous, and the result is
$$
y(x)\sim e^{-x^2/2}x^{\nu}\,,
$$
where $\nu=\epsilon-1/2$. At this point, we guess that the solution can be written as
$$
y(x) = e^{-x^2/2}x^{\nu}\sum_{n=0}^{\infty}a_nx^{-n}\,.
$$
(Note that this is a power series in negative powers of $x$!) This results in a recurrence relation for the $a_n$'s, given by
$$
a_1=0,~~~~~
a_{n+2} = -\frac{(n-\nu)(n-\nu+1)}{2(n+2)}a_n\,,
$$
so that the solutions looks like
$$
y(x) = e^{-x^2/2}x^{\nu}
\left(
1-\frac{\nu(\nu-1)}{4x^2} + \frac{\nu(\nu-1)(\nu-2)(\nu-3)}{32x^4}
- \frac{\nu(\nu-1)(\nu-2)(\nu-3)(\nu-4)(\nu-5)}{384x^6} + \cdots
\right)
$$
Crucially, we can compute the radius of convergence of this series, and it is zero. It is still a good asymptotic solution for any $\nu$, but it is not good for the situation we are concerned with, which is for the function to go to zero at infinity. This computation of the radius of convergence of the asymptotic series replaces the hand-wavy "the power
series acts like $e^{x^2}$" with a slightly less hand-wavy "the radius of convergence of the series expansion about $x\to\infty$ is zero".
In any case, this means that we only get a "good" solution when the power series truncates to a finite number of terms (because in that case the radius of convergence is of course infinite!), and this occurs when $\nu$ is a non-negative integer, which can be clearly seen from both the recursion relation and the first few terms of the expansion. After multiplying through by the factor of $x^{\nu}$, we get exactly the physicist Hermite polynomials $H_n(x)$, as we should.
The problem is that it's not clear to what extent the radius of convergence being zero is really related to the function blowing up at infinity, so this doesn't entirely answer the question, unfortunately.
Best Answer
This is for $k=1$, for which your formulas are valid.
In short, they come from Fourier integral representations of solutions to the ODE. This is much easier to see for $Ai$ than $Bi$.
Applying the standard Fourier transform $$ \tilde\Psi(\omega) = \mathcal F\{\Psi\} = \int_{-\infty}^\infty \Psi(x)e^{-i\omega x}dx $$ to the ODE, leads to the first order ODE $-\omega^2\tilde\Psi - i \frac{d\tilde\Psi}{d\omega} = 0$, thus $$ \tilde\Psi = Ae^{i\omega^3/3}. $$ Substituting into the inverse transform results in the above form $$ \Psi(x) = \mathcal F^{-1}\tilde\Psi = \frac{A}{2\pi}\int_{-\infty}^\infty e^{i\omega^3/3 + ix\omega} d\omega $$ (the imaginary part cancels over the entire line leading to your formula, while the constant $A$ is chosen with respect to the far field behaviour).
This has only recoverd $Ai$, rather than $Bi$, since implicit in applying the Fourier transform is the idea that the solution is integrable on $(-\infty,\infty)$, which $Bi$ is not.
The following approach to find other solutions is due to Ablowitz and Fokas (Complex Variables: Introduction and Applications, 2nd ed). Consider a more general Fourier-like integral representation of a solution to Airy's equation $$ \Psi(x) = \int_C V(\omega)e^{i\omega x} d\omega $$ where $C$ is a complex contour whose properties are yet to be decided (but we think of going to infinity in both direction along certain rays). Using integration by parts we have $$ \Psi'' - x\Psi = \int_C [-\omega^2 V - i V'(\omega)]e^{i\omega x} d\omega - i[V(\omega)e^{i\omega x}]_C. $$ Here the latter term represents the contribution from each endpoint of $C$, and the integration by parts is valid since $V$ will end up being entire. To set this to zero we have, as before, $V = i\omega^3/3$, and need to take $C$ such that $[e^{i(\omega^3/3 + \omega x)}]_C = 0$. If $\arg\omega=\theta$, this works for $\omega \to \infty$ in the sectors $$ \sin(3\theta) > 0 $$ or $0 < \theta < \pi/3$, $2\pi/3 < \theta < \pi$, $-2\pi/3 < \theta < -\pi/3$.
Now the Fourier integral solution $Ai$ can be thought of as the result when $C$ is the curve connecting $\theta = \pi$ to $\theta = 0$ (in a limiting sense, as these angles are on the border between exponential decay and growth, hence the oscillatory nature of the integrand). However, if one takes $C$ to connect $\theta = -\pi/2$ to $\theta = 0$, that is, $C$ comes up the imaginary axis from $-i\infty$ and then goes out along the real axis to $+\infty$, the integral ends up being (taking $\omega = -i\zeta$ on $\theta = -\pi/2$) \begin{align*} \Psi(x) &= i\int_0^\infty e^{i(i\zeta^3/3 - ix\zeta)} d\zeta + \int_0^\infty e^{i(\omega^3/3 + x\omega)} d\omega \\ &= \int_0^\infty \cos(\omega^3/3 + x\omega) d\omega + i\int_0^\infty \sin(\omega^3/3 + x\omega)d\omega + i\int_0^\infty e^{-\zeta^3/3 + x\zeta} d\zeta. \end{align*} This result is of course up to the arbitrary multiplicative constant in $V$. The first term is proportional to $Ai$, so if it is subtracted off and the arbitrary constant is chosen appropriately (again with respect to some far field behaviour), you end up with your integral formula for $Bi$.