following the discussion with the OP:
Basically you are right - we just need a transformation from one basis to another. Namely, from the eigenstates of $H$ before the turning on of the field to the eigenstates of $H$ after. Let's call them $H$ and $H'$ for convenience, and their respective eigenstates $|n\rangle$ and $|n'\rangle$. As you notice, this transformation is achieved by a unitary matrix $U$, that relates the expansion coefficients in one basis $|\psi \rangle = \sum_n a_n |n\rangle$ to the expansion coefficients in another basis $|\psi \rangle = \sum_{n'} a_{n'} |n'\rangle$. So $a_{n} = \sum_n U_{n, n'} a_{n'}$. 66
It is easy to see that the coefficients are the overlap between the basis states, as $U_{n, n'} = \langle n | n' \rangle$. To see this, just apply it to a basis set $|n'\rangle = \sum_{n} U_{n, n'} |n\rangle$ and enter the unit operator $\sum_{m}|m\rangle \langle m|$ to get $$ \sum_{m}|m\rangle \langle m|n'\rangle = \sum_{n}U_{n,n'}|n\rangle $$ which from the orthogonality of the basis state will give you $U_{n,n'} = \langle n | n'\rangle$. This overlap is exactly the integral over the two wave functions.
However, the integral is not the only way to calculate the overlap. As you noted, the basis of $H'$ can be obtained from the basis of $H$ by a translation in $x$. This translation can be written as an operator $$T_{x_0} = e^{-ipx_0/\hbar} = e^{\sqrt{\frac{m \omega}{2\hbar }}x_0(a-a^{\dagger})} = e^{-\frac{m \omega x_0^2}{4\hbar }}e^{-\sqrt{\frac{m \omega}{2\hbar}}x_0a^{\dagger}}e^{\sqrt{\frac{m \omega}{2\hbar }}x_0a}$$ so for the ground state the overlap is given by $$\langle 0' | 0 \rangle = e^{-\frac{m \omega x_0^2}{4\hbar }} \langle 0 |e^{-\sqrt{\frac{m \omega}{2\hbar }}x_0a^{\dagger}}e^{\sqrt{\frac{m \omega}{2\hbar }}x_0a} |0 \rangle = e^{-\frac{\hbar m \omega x_0^2}{4}}$$ where we used the fact that $a|0\rangle = 0$ so only the zeroth term in the power expansion of the exponent can contribute 66