Why do we only look at transitions between energy-eigenstates, when perturbing a system with an oscillating interaction? I have some doubts about the usual explanation of sharply peaked emission / absorption spectra, that one can observe when one looks at quantum mechanical systems, for example the hydrogen atom.
In short: We can calculate the transition probability between two energy eigenstates, and we see that for a given time, its amplitude is sharply peaked exactly at the energy difference between the states. That is clear. But why are "transitions" between energy eigenstates the only ones that would occur in nature in general?
In more detail:
I look at a system in the Schrödinger-picture where the time evolution is given by $H_0$. Now I perturb the system, using an oscillating time dependent interaction $H_{\text{int}} = \lambda e^{-i\omega t} H_i$.
For little $\lambda$, we can expand the solution in $\lambda$ (Dyson Series) and calculate the probability-amplitude that the system, prepared as $|n\rangle$ at $t_0$, can be found as $|m\rangle$ at time $t$. To measure this, we look at the expectation value of the projection operator $P = |n \rangle \langle n |$. Then
\begin{align}
\langle P \rangle_{|n(t)\rangle} = | \langle m | U_{\omega}(t, t_0) | n \rangle |^2
\end{align}
With $U_{\omega}$ being the solution to
\begin{align}
i \hbar \partial_t U_{\omega}(t, t_0) = (H_0 + H_{\text{int}, \omega} )U(t, t_0)
\end{align}
Employing the Dyson-series, we can calculate this probability for a given time t, and a given frequency $\omega$ of the interaction, to an order $\lambda^n$ in the strength of the interaction. The result will be sharply peaked for $\omega$ being close to $\omega_m - \omega_n$. For the interaction being subject to a time-dependent electric field (e.g. light), we interpret a change in the states as something that changes the electric field, and thus some effect that can be measured while changing the fields electricity.
As such, we can measure those sharp peaks in an absorption or emission spectrum, of atoms for example.
Now - Why can't we do the same calculation for abitrary states? Nothing prevents me from writing down, and calculating:
\begin{align}
| \langle \Psi | U_{\omega}(t, t_0) | \Phi \rangle |^2
\end{align}
For abitrary $ | \Psi \rangle $, and $|\Psi \rangle $, and this will be a correct result. But I won't see abitrary peaks in my emission / absorption spectrum, but only the peaks from the energies.
Is it because I can expand
\begin{align}
| \Psi \rangle = \sum_m \Psi_m(t) | m \rangle
\end{align}
and
\begin{align}
| \Phi \rangle = \sum_n \Phi_n(t) | n \rangle
\end{align}
and thus
\begin{align}
 \langle \Psi | U_{\omega}(t, t_0) | \Phi \rangle |^2 = | \sum_{m, n} \Psi_m(t)^* \Phi_n(t) \langle m | U_{\omega}(t, t_0) | n \rangle|^2
\end{align}
I however don't see why the last expression should still be sharply peaked around distinct $\omega$-values. It's a double countable infinite sum over sharply peaked functions of $\omega$, with possibly time dependent coefficients.
Am I on the right track here, or is there another way to see this?
 A: You can't "see" the psi function (or expansion coefficients) having strong response to frequency of external wave $\Omega$ approaching difference frequencies $\omega_{nk} = \frac{E_n - E_k}{\hbar}$ from the general formulae above, as you have demonstrated. They are too general to a point where no assumption means no conclusion.
You really have to make some assumptions, such as starting with some Hamiltonian eigenfunction that oscillates sinusoidally without a perturbation, and calculate response due to external oscillating electric field (and then maybe see how this response depends on external frequency $\Omega$).
Erwin Schroedinger did this approximately in his paper Quantisation as a Problem of Proper Values (Part IV), Annalen der Physik (4), vol. 81, 1926. See the English translation in the book Collected Papers on Wave Mechanics by E. Schroedinger, Blackie & Son Ltd., 1928.
He found out that (system starting in $k$-th eigenfunction of $H_0$), the resulting reaction of $\psi$ to oscillating perturbation is approximately the combination of two oscillations with frequencies $E_k/\hbar - \Omega, E_k/\hbar + \Omega$, and with amplitude (spatial part) that is proportional to (assuming $\Omega$ does not equal any of the difference frequencies):
$$
w_{\pm}(\mathbf r) = \sum_{n=1}^\infty \frac{H'_{nk} \Phi_n(\mathbf r)}{E_n - E_k \pm \hbar \Omega}.
$$
where $H_{nk}'$ is matrix element of the perturbation Hamiltonian and $\Phi_n$'s are eigenfunctions of the unperturbed Hamiltonian $H_0$.
This means that when external frequency $\Omega$ gets close to some difference frequency $\omega_{kn}$ implied by $H_0$, correction to $\psi$ gets much greater (due to a very small number in the denominator of the implied term of the series).
One could "guess" there will be this kind of resonant response to external oscillation even without carrying out the perturbation theory calculations. For non-eigenfunction $\psi$, oscillation of dipole moment during "free oscillation", governed by $H_0$, is sum of oscillations at difference frequencies $\omega_{kn}$ with each oscillation being as strong as expansion coefficients and dipole moment matrix elements (in $H_0$ basis) determine. Exchange of energy between field and matter can happen systematically only when electric current oscillation (due to dipole oscillating) is in phase with oscillation of external electric field, and this requires very similar frequencies of both oscillations, otherwise in the long run of many periods, there is no systematic transfer and energy just oscillates between matter and field. So as long as the free oscillation is not destroyed by the perturbation too much, systematic absorption requires that external wave have frequency close to some of the $\omega_{nk}$'s. However, this is just a guess, because without doing the calculations, we don't have a good estimate of what the free oscillation changes into under the perturbation.
Schroedinger also derives from his calculations expected average of electric dipole moment in direction of axis $y$ (surprisingly, he chose a direction perpendicular to that of electric field):
$$
\langle \boldsymbol {\mu}_y \rangle = a_{kk} + 2E_0\sum_{n=1}^{\infty} \frac{(E_n - E_k)a_{kn}b_{kn}}{(E_n - E_k)^2 - \hbar^2 \Omega^2} \cos \Omega t
$$
where $a_{kn}, b_{kn}$ are some matrix elements and $E_0$ is amplitude of external electric field. Similar formula would probably come out for the $x$ component, relevant for absorption.
Thus dipole moment is forced to oscillate at the external frequency $\Omega$ (as it would in a simple linear classical model), but in addition, amplitude of its oscillations depends on external frequency as if there are multiple resonant frequencies.
The model is deficient in that there is no way for the absorbed energy to dissipate and for the molecule to reach some kind of steady state in the external field. A more detailed calculation would probably show this model goes more and more into higher energy states, something that does not seem realistic if external frequency is below the ionization limit. In reality, spontaneous emission and other interactions with environment would keep the higher states unoccupied, but this is not taken into account here.
For a toy model of a two-level system, one can probably do calculations of these things without approximations of the perturbation theory, instead doing a full numerical solution. See the Rabi oscillations.
A: 
Why do we only look at transitions between energy-eigenstates, when perturbing a system with an oscillating interaction?

We don't have to.
When we calculate transition rates we typically sum over all final states:
$$
R \sim \sum_{f}|\langle \Psi_f|V|\Psi_0\rangle|^2\delta(E_0 + \omega - E_f)\;.
$$
The sum can be performed in any basis. But, why make your life difficult?
The initial state $\Psi_0$ can be any state, but how are you going to prepare a state other than a stationary state (energy eigenstate)? If the state is not stationary, how are you going to make it last long enough to hit it with a beam of probe particles? In many cases the state $\Psi_0$ has to be the ground state because there is no other way to prepare a system that won't decay before you hit it with the beam.
