Broadband light term in a Hamiltonian In atomic systems, for a two-level system, the Hamiltonian can be written in the form:
$$H=\left(
\begin{array}{cc}
 E_1 & C_{12} \\
 C_{21} & E_2 \\
\end{array}
\right)$$
where $E_1$ and $E_2$ are the energy offsets of the atomic levels from the zero point, and $C_{**}$ are the coherences between the levels. Now if we would like to couple the two levels with an electromagnetic field, usually the coupling term looks like this:
$$C_{12}=\Omega \cos \omega t$$
where $\Omega$ is the Rabi frequency of the system and $\omega$ is the frequency of the electric field, or the Larmor frequency.
The question: Now this kind of light-coupling couples light with zero linewidth. No broadband light. So my question is, how can we include light that has a linewidth? 
More details:
Now normally I include states that have a width not using the Hamiltonian, but using the Liouville equation. It looks like this:
$$i\hbar\frac{d\rho}{dt}=H \rho-\rho H-\frac{1}{2}i \hbar(\Gamma\rho+\rho\Gamma)$$
where $\rho$ is the density matrix of the system. $\Gamma$ is a diagonal matrix that contains the linewidth of each state. Is there a similar solution to include the width of the electromagnetic field? That would be OK too.
Why do I need this? I have a complicated, large system with many levels and I would like them all to be excited together with a broadband light, kind of like the experiment I'm working on.
 A: That attempt to account for dissipation in the von Neumann equation hasn't been in use for about 40 years now: it doesn't work very well. First and foremost, it "leaks" probability, since $Tr\dot\rho \neq 0$. Second, but probably most importantly, it cannot keep the density matrix positive definite. 
If you want a dissipative evolution in terms of the two-level density matrix only, then what you are looking for is a (linear) Lindblad dynamics of the form
$$
\dot\rho = -i[H,\rho] + \frac{1}{2}\sum_j{\left(2\Gamma_j\rho\Gamma_j^\dagger - \Gamma_j^\dagger\Gamma_j\rho - \rho\Gamma_j^\dagger\Gamma_j\right)}
$$ 
In this case, the density matrix $\rho$, the Hamiltonian $H$, and the dissipative operators $\Gamma_j$ all act on the two-level Hilbert space only. The Hamiltonian $H$ and the $\Gamma_j$-s can also be time-dependent. This is the most general form of a positive linear dynamics, i.e. one that keeps $\rho$ positive definite at all times. You can easily check that it conserves total probability too, in the sense that $Tr\dot\rho = 0$. 
A frequent choice for the dissipative operators is $\Gamma_j = \sigma_\pm$, but in the simplest possible case, the sum over $j$ keeps a single term and the evolution eq. is simply
$$
\dot\rho = -i[H,\rho] + \frac{1}{2}\left(2\Gamma\rho\Gamma^\dagger - \Gamma^\dagger\Gamma\rho - \rho\Gamma^\dagger\Gamma\right)
$$
This can actually be solved exactly. If you rewrite your $\Gamma$ as $\Gamma \rightarrow \Gamma^\dagger\Gamma$, then the form above differs from yours only by a $\Gamma\rho\Gamma^\dagger$ term.  
If you prefer to account explicitly for the electromagnetic field, it can be accommodated directly in the Lindblad dynamics, provided the Hilbert space is extended to the system-field Hilbert space, and the density matrix becomes the system-field density matrix. In this case the total hamiltonian must include the two-level hamiltonian $H_0$, the field hamiltonian $H_1$, and the system-field interactions,  $H_{int}$, so 
$$
H = H_0 + H_1 + H_{int}
$$
The dissipative terms usually account for possible field interactions with another external bath/reservoir and depend only on the field degrees of freedom. 
The system-field hamiltonian is known in general as the Jaynes–Cummings model and can be found in any course on Quantum Optics. I am jotting it down here just to make the answer self-contained. First write the two-level hamiltonian in terms of SU(2) operators,
$$
H_0 = \epsilon \hat\sigma_z + V\hat\sigma_+ + V^*\hat\sigma_-
$$
where $\hat\sigma_z = |e\rangle\langle e| -  |g\rangle\langle g|$, $\hat\sigma_+ = |e\rangle\langle g|$, and $\hat\sigma_- = |g\rangle\langle e|$, with $|e\rangle$ the excited state and $|g\rangle$ the ground state. If you prefer, it is also possible to write instead
$$
H_0 = \vec{u} \cdot\vec{\hat\sigma}
$$
where $\vec{\hat\sigma} = (\hat\sigma_x, \hat\sigma_y, \hat\sigma_z)$ and $\vec{u}$ is a real 3-D vector (possibly a magnetic field?). 
The field hamiltonian is just the standard hamiltonian for arbitrary field modes,
$$
H_1 = \sum_{\vec{k},\lambda}{\omega(\vec{k})a^\dagger_{\vec{k},\lambda}a_{\vec{k},\lambda}}
$$
and the system field interaction is initially taken as (Schroedinger representation)
$$
H_{int} = \sum_{\vec{k},\lambda}{(g_{\vec{k},\lambda}\sigma_+ + g^*_{\vec{k},\lambda}\sigma_-)(a^\dagger_{\vec{k},\lambda} + a_{\vec{k},\lambda})}
$$
but after a rotating wave approximation (RWA) in the interaction picture and reverting back to Schroedinger it becomes
$$
H_{int} = \sum_{\vec{k},\lambda}{(g_{\vec{k},\lambda}\sigma_+a_{\vec{k},\lambda} + g^*_{\vec{k},\lambda}\sigma_-a^\dagger_{\vec{k},\lambda})}
$$
As for the dissipative operators $\Gamma_j$ they are usually defined as 
$$
\Gamma_j \;\;\rightarrow \;\; \Gamma_{\vec{k},\lambda} = \kappa_{\vec{k},\lambda}a_{\vec{k},\lambda} 
$$
where $\kappa_{\vec{k},\lambda}$ are complex scalars (coupling constants).
Finally, the broadband condition is not so much about the form of the evolution equation as about the initial state of the field. If the field is initially confined to a single mode, the system is basically interacting with a coherent laser mode. If instead the field starts in some pure multi-mode superposition, or even better, some multi-mode mixed state, then you basically end up with a broadband field.
