Edited after comments.
I can't seem to find (or derive) the exact solution of the two quantum state system absorption problem in the following formulation. Let's consider the general harmonic wave absorption:
$$i \hbar \partial_t \Psi(\mathbf{r},t)=\left(\hat{H_0}(\mathbf{r})+\hat{H}(\mathbf{r}) e^{-\gamma |t|} \cos \omega t \right)\Psi(\mathbf{r},t) \tag{1}$$
Where the optical pulse has the exponential decay both for positive and negative time (it's the model I decided to use for simplicity).
Using the eigenstates of $\hat{H_0}$, we expand the wavefunction:
$$\Psi(\mathbf{r},t)=\sum_n C_n(t) \psi_n(\mathbf{r}) e^{-i \omega_n t} \tag{2}$$
Using orthonormality ($H_{mn}$ are matrix elements of the harmonic wave amplitude $\hat{H}(\mathbf{r})$):
$$i \hbar \dot{C_m}(t)=\sum_n C_n(t) H_{mn} \cos \omega t~ e^{i (\omega_m-\omega_n) t-\gamma |t|} \tag{3}$$
Now assume that we have a simple two state system.
$$\begin{bmatrix} \dot{C_1} \\ \dot{C_2} \end{bmatrix}=\frac{e^{-\gamma | t|} \cos \omega t}{i\hbar}\begin{bmatrix} H_{11} & H_{12} e^{i \omega_{12} t} \\ H_{12} e^{-i \omega_{12} t} & H_{22} \end{bmatrix} \begin{bmatrix} C_1 \\ C_2 \end{bmatrix} \tag{4}$$
How to obtain the exact general solution for (4)? A reference would be enough for me, because I was only able to find either some particular case or a perturbation treatment of the problem.
Here's a related question but with a different matrix.
What's interesting, if we use Fourier transform:
$$C_n=\frac{1}{\sqrt{2 \pi}} \int_{-\infty}^\infty S_n(u) e^{-i u t} du \tag{5}$$
Then we obtain a pair of coupled integral equations:
$$\frac{2 \pi}{\gamma}\hbar u S_1(u)=H_{11}\int_{-\infty}^\infty S_1(v) \left( \frac{dv}{(v-u-\omega)^2+\gamma^2}+\frac{dv}{(v-u+\omega)^2+\gamma^2} \right)+ \\ +H_{12}\int_{-\infty}^\infty S_2(v) \left( \frac{dv}{(v-u-\omega+\omega_{12})^2+\gamma^2}+\frac{dv}{(v-u+\omega+\omega_{12})^2+\gamma^2} \right) \tag{6}$$
$$ \frac{2 \pi}{\gamma}\hbar u S_2(u)=H_{22}\int_{-\infty}^\infty S_2(v) \left( \frac{dv}{(v-u-\omega)^2+\gamma^2}+\frac{dv}{(v-u+\omega)^2+\gamma^2} \right) + \\ + H_{12}\int_{-\infty}^\infty S_1(v) \left( \frac{dv}{(v-u-\omega-\omega_{12})^2+\gamma^2}+\frac{dv}{(v-u+\omega-\omega_{12})^2+\gamma^2} \right) \tag{7}$$
Where the initial conditions for $C_n$ uniquely determine the normalization, as could be seen from (5).
These should be easy enough to solve numerically and then use FFT to get back to $C_n(t)$.
If this problem was never treated exactly or numerically, then I suppose it could be fun to try it on my own and see how the results are different from the usual approximations.
