I am interested in a rigorous derivation of the coherent evolution of a quantum system with $N$ states under the application of a harmonic perturbation. I have read several articles about the use of adiabatic elimination for this, but none of them seem to be what a mathematician would call rigorous. Some just set a time derivative to zero, some assume constant amplitudes, some make a rotating wave approximation even before adiabatically eliminating some states.
Assume a state vector $|\psi(t)\rangle$ in an n-dimensional Hilbert space that evolves under $$i\hbar\frac{d}{dt}|\psi(t)\rangle_S = \left(\hat{H}_0 + \hat{V}_S\cos(\omega_p t + \phi) \right) |\psi(t)\rangle_S$$ where $H_0$ is a time-independent Hamiltonian and $\omega_p$ is the frequency of the applied perturbation. This is often due to interaction of a quantum system with radiation.
In the interaction picture, defined by the transformation $|\psi(t)\rangle_I = e^{i\frac{\hat{H}_0}{\hbar} t}|\psi(t)\rangle_S$, the evolution is then given by $$i\hbar\frac{d}{dt}|\psi(t)\rangle_I = \hat{V}_I\cos(\omega_p t + \phi)|\psi(t)\rangle_I$$ where $\hat{V}_I = e^{i\frac{\hat{H}_0}{\hbar} t} \, \hat{V}_S \, e^{-i\frac{\hat{H}_0}{\hbar} t}$.
Projecting this onto the basis $\{|i\rangle\}$ of eigenstates of $\hat{H}_0$, with ${H}_0|i\rangle = E_i |i\rangle$, this leads to coupled equations in the interaction picture: $$\frac{d c_i(t)}{dt} = -\frac{i}{2} \sum_j \; \Omega_{ij} \; c_j(t) \; \left(e^{i[(\omega_{ij} + \omega_p)t + \phi]} + e^{i[(\omega_{ij} - \omega_p)t - \phi]}\right)$$ where I have defined $\Omega_{ij} = \frac{\langle i|V_S|j\rangle}{\hbar}$ (which is a measure of the strength of the perturbation connecting two states $|i\rangle$ and $|j\rangle$) and $\omega_{ij} = \frac{E_i - E_j}{\hbar}$.
How can one show with mathematical rigour that under certain constraints for $\Omega_{ij}$ and $\omega_p$ relative to $\omega_{ij}$, in combination with initial values for $c_i(t=0)$, only certain states will have a slow evolution and will be coupled effectively? In particular, if $\omega_p$ is close to exactly one $\omega_{ij}$, I expect that the system can effectively be reduced to a two-level system, with a shift in energy (for electric fields, this is the AC Stark shift) by the effect of off-resonant levels. The fact that the population of off-resonant states will be small should come out of such an adiabatic elimination (integrating out the fast variables), rather than be put in at the beginning.
In other words, what is a rigorous way of deriving an effective Hamiltonian involving a subset of states from the general Schrödinger equation as given above?
