To what extent can one recover plane waves from the Airy eigenfunctions of a linear potential as the field is turned off? Consider a single massive particle in one dimension under the action of a static linear potential, with the hamiltonian
$$
\hat H=\frac{\hat p^2}{2}+\hat{x}F_0.
$$
The eigenstate at energy $E$ is, with this normalization, given by
$$
\langle x|\chi_E\rangle
=
\frac{2^{1/3}}{F_0^{1/6}}
\operatorname{Ai}\left(\sqrt[3]{2F_0}\left(x-\frac{E}{F_0}\right)\right),
$$
where $\operatorname{Ai}$ is the Airy function and the eigenstates satisfy $\int_{-\infty}^\infty|\chi_E\rangle \langle \chi_E|\mathrm dE=1$ and $\langle \chi_E|\chi_{E'}\rangle=\delta(E-E')$.
I would like to know to what extent I can switch off the field, i.e. investigate the limit $F_0\to 0$, and recover the plane-wave eigenstates of the free particle.
Much of this is routine: I can apply the Airy function's asymptotic property, 
$$
\operatorname{Ai}(-z)\sim\frac{1}{\sqrt{\pi}z^{1/4}}\sin\left(\tfrac23 z^{3/2}+\frac \pi 4\right),
$$
plus some elementary Taylor series on the resulting powers of $1-F_0 x/E$ to get that, if $E>0$ is fixed and $x$ is in a fixed, bounded interval, then
$$
\langle x|\chi_E\rangle
\sim
\frac{2^{1/4}}{\sqrt{\pi}E^{1/4}}
\sin\left(
\frac\pi4 +\frac{\sqrt{8}}{3}\frac{E^{3/2}}{F_0}-\sqrt{2E}x
\right)
. \tag{1}
$$
This is, all things told, pretty good. I recover the plane waves with the appropriate momentum $|p|=\sqrt{2E}$, and they're even pretty much correctly normalized (since $\mathrm dE/(\sqrt[4]{E})^2\propto \mathrm dp$). Since I'm taking the limit of real-valued functions, I naturally get a real-valued limit.
However, this isn't perfect. For one, I only recover one linearly independent eigenstate at each energy. There should be another, phase-shifted solution, which would be proportional to a cosine of the same argument. For a very small but nonzero $F_0$, this cosine solution will turn into an Airy function of the second kind $\operatorname{Bi}$ at about $x=E/F_0$, and then blow up super-exponentially after that. However, my analysis in a pre-chosen bounded interval of $x$ and sufficiently small $F_0$ can't know anything about what will happen way over there, nor should it really care.
Most importantly, though, is the fact that the limit is not really well defined, because of that horrible extra phase in $E^{3/2}/F_0$, which makes everything not tend to anything. You could say that in a way this solves my previous problem, because you can write the approximation as
\begin{align}
\langle x|\chi_E\rangle
\sim
\frac{2^{1/4}}{\sqrt{\pi}E^{1/4}}
&
\left[
\sin\left(
\frac\pi4 +\frac{\sqrt{8}}{3}\frac{E^{3/2}}{F_0}
\right) \cos\left(
\sqrt{2E}x
\right)
\right. \\ & \qquad \left.
-
\cos\left(
\frac\pi4 +\frac{\sqrt{8}}{3}\frac{E^{3/2}}{F_0}\right) \sin\left(\sqrt{2E}x
\right)
\right].
\end{align}
If I fix $E$ and a small but nonzero $F_0$ such that the solution behaves like the first, cosine term, then there is another eigenstate at energy $E'=E+\delta E$ slightly above it such that
$$
\frac{\sqrt{8}}{3}\frac{E'^{3/2}}{F_0}=\frac{\sqrt{8}}{3}\frac{E^{3/2}}{F_0}+\frac\pi2
$$
or, to first order in $F_0$, $\delta E=\frac{\pi}{2}\frac{F_0}{\sqrt{2E}}$. As $F_0$ goes to zero $\delta E$ also shrinks and this second solution comes to sit at the same energy as I started with, but with the nontrivial phase that I needed.
This analysis mostly strikes me as disingenuous, and it definitely ignores the fact that the wave at fixed energy in $(1)$ does not actually converge to anything and it does not have a well-defined phase if it does. This second, mysterious solution at some other energy $E'$, which depends on $F_0$, will also not converge to anything nor will it have a definite phase, but this non-definite phase will somehow mysteriously be exactly complementary to the non-definite phase I started with. Phrased in that light, it makes no sense at all - though to be honest it does feel like it has enough of the ingredients of nontriviality that it can be built up into an argument that's actually coherent.

So: I am looking for methods or references to deal rigorously with the $F_0\to0$ limit, and to extract by rigorous means, if possible, the full set of two linearly independent eigenstates per definite energy of the free particle from this system.
 A: This is meant as a comment but too long. Hence I post it as an
answer.
One can also look at the time evolution and consider the $F_{0}\rightarrow 0$
limit in that context. The commutation properties of the operators involved
allow us to rewrite the time evolution operator as a product where each
factor no longer contains non-commuting operators. Thus we have
\begin{equation*}
\exp [-iHt]=\exp [-i(\frac{p^{2}}{2}+F_{0}x)t]
\end{equation*}
and find by differentiation wrt time that
\begin{eqnarray*}
\partial _{t}U(t) &=&\partial _{t}\{\exp [-iF_{0}xt]\exp [-iHt]\}=\exp
[-iF_{0}xt]\{-i(F_{0}x-\frac{p^{2}}{2}-F_{0}x)\}\exp [-iHt] \\
&=&=\frac{-i}{2}\exp [-iF_{0}xt]p^{2}\exp [+iF_{0}xt]U(t) \\
&=&=\frac{-i}{2}(p+F_{0}t)^{2}U(t)
\end{eqnarray*}
so
\begin{eqnarray*}
U(t) &=&\exp [\frac{-i}{2}(p^{2}t+pF_{0}t^{2}+\frac{1}{3}F_{0}^{2}t^{3})] \\
\exp [-iHt] &=&\exp [+iF_{0}xt]\exp [\frac{-i}{2}(p^{2}t+pF_{0}t^{2}+\frac{1
}{3}F_{0}^{2}t^{3})]
\end{eqnarray*}
The advantage is that the exponentials no longer contain non-commuting
operators . We require strong convergence
\begin{equation*}
\exp [-iH(F_{0})t]f\overset{F_{0\rightarrow 0}}{\rightarrow }\exp
[-iH(0)t]f=\exp [\frac{-i}{2}p^{2}t]f
\end{equation*}
We can try to estimate
\begin{eqnarray*}
&\parallel &\exp [-iH(F_{0})t]f-\exp [-iH(0)t]f\parallel \\
&=&\parallel \exp [+iF_{0}xt]\{\exp [\frac{-i}{2}(p^{2}t+pF_{0}t^{2}+\frac{1%
}{3}F_{0}^{2}t^{3})]f-\exp [\frac{-i}{2}(p^{2}t]\}f\parallel \\
+ &\parallel &\{\exp [+iF_{0}xt]-1\}\exp [\frac{-i}{2}(p^{2}t]f\parallel \\
&\leqslant &\parallel \{\exp [\frac{-i}{2}(p^{2}t+pF_{0}t^{2}+\frac{1}{3}%
F_{0}^{2}t^{3})]f-\exp [\frac{-i}{2}(p^{2}t]\}f\parallel \\
+ &\parallel &\{\exp [+iF_{0}xt]-1\}\exp [\frac{-i}{2}(p^{2}t]f\parallel \\
&=&\parallel \{\exp [\frac{-i}{2}(pF_{0}t^{2}+\frac{1}{3}F_{0}^{2}t^{3})]f-f%
\parallel +\parallel \{\exp [+iF_{0}xt]-1\}\exp [\frac{-i}{2}%
(p^{2}t]f\parallel \rightarrow 0
\end{eqnarray*}
This seems all right, i.e. the limit makes sense. However it does not tell
us anything about the Airy functions. Finding corrections to leading order
in $F_{0}$ looks challenging but quite difficult.
A: Regarding the first part of your question. It is a matter of order of limits. In your setup you have (at least) two parameters: $F_0$ and $x$. Saying that $F_0 \to 0$ is not enough, you have to specify what happens with say $F_0 x$. Clearly, for whatever small (but nonzero) $F_0$ if you allow $x$ to be arbitrary, then the product $F_0 x$ can be arbitrary large. In this case there is no reason to expect that the initial problem be reducible to a free particle one, for the perturbation $F_0 x$ is large. If you want to recover plane waves you first consider the limit $F_0 \to 0$ for other parameters fixed, and later you may study $x\to \pm \infty$ as well.
The same goes for the second part: both $Ai(z)$ and $Bi(z)$ are oscillating functions when $z\to -\infty$, which is precisely the limit you need. You have all the ingredients. Two solutions in the limit $F_0\to 0$ behave as
\begin{eqnarray}
\psi^E_1(x) & = & \frac{A_1}{E^{1/4}}\sin\left (\frac{c}{f}+\frac{\pi}{4}-\sqrt{2E}x \right ) \\
\psi^E_2(x) & = & \frac{A_2}{E^{1/4}}\cos\left (\frac{c}{f}+\frac{\pi}{4}-\sqrt{2E}x \right )
\end{eqnarray}
Instead of dealing with sines and cosines we can form a linear compination
\begin{eqnarray}
\tilde \psi^E_1(x) & = & \frac{\tilde A_1}{E^{1/4}}\exp\left (\frac{c}{f}+\frac{\pi}{4}-\sqrt{2E}x \right ) = \frac{B_1}{E^{1/4}}\exp\left (-\sqrt{2E}x \right ) \\
\tilde\psi^E_2(x) & = & \frac{\tilde A_2}{E^{1/4}}\exp\left (-\frac{c}{f}-\frac{\pi}{4}+\sqrt{2E}x \right )=\frac{B_2}{E^{1/4}}\exp\left (\sqrt{2E}x \right ),
\end{eqnarray}
where we just absorbed an unimportant "divergent phase" in the normalization constant
A: 
Most importantly, though, is the fact that the limit is not really well defined, because of that horrible extra phase in $E^{3/2}/F_0$, which makes everything not tend to anything.

This ill-definedness should be expected: as you reduce $F_0$, the classical turning point goes farther from $x=0$ (assuming $E>0$), continuously. The phase will of course also continuously run away.
So, the most you can do is extract this phase and take the limit for the remaining part of the wavefunction. Since you accept the use of the divergent $\operatorname{Bi}$ function, why not just start from the running wave $\operatorname{Ai}(x)+i\operatorname{Bi}(x)$ instead of the separate standing waves? In this case you'll get the asymptotics
$$\frac{2^{1/4}}{\sqrt{\pi}E^{1/4}}\exp\left(i\left(-\frac\pi4 +\frac{\sqrt8}3\frac{E^{3/2}}{F_0}-\sqrt{2E}x\right)\right),$$
from which the running $F_0$-dependent phase can be conveniently extracted as a global factor. After doing this you're left with the desired plane wave.
