Two-level system with spontaneous emission I have been stuck with this elementary problem of two-level system including spontaneous decay. After solving by standard procedure the following pair of equations, when I plot the expressions of the ground state population $P_0 (=|C_0(t)|^2)$, excited state population $P_1 (=|C_1(t)|^2)$ and sum of the two, $P_0+P_1$ vs time $t$, I find that $P_0+P_1$ is not $1$, as if there is a loss of the system as a whole. It has got to be wrong.
I tried to solve the following coupled differential equations incorporating the spontaneous decay width, for a single two-level atom subjected to optical field:
$ \dot{C_0(t)} = -\frac{i}{2} \Delta C_0(t) + \frac{i}{2}(\Omega-i\Gamma_{sp})C_1(t) \\
 \dot{C_1(t)} = \frac{i}{2} (\Delta+i\Gamma_{sp}) C_1(t) + \frac{i}{2}(\Omega^*)C_0(t) $
Is it correct to incorporate spontaneous emission in the language of probability amplitude of a single atom like this?
 A: This is not really the way to account for spontaneous emission in a two-level system. The reason for this is that emitting a photon will land you in an entangled state with the field, which means that the coherence of the atomic part of the system is diminished, and that part cannot be considered to be in a pure state. To to this, then, you need to describe your system in terms of a density matrix rather than a pure state, and the resulting formalism is known as the optical Bloch equations.
To summarize them here, what you do is work in terms of the matrix
$$
\rho
=\begin{pmatrix}\rho_{gg}&\rho_{ge}\\ \rho_{eg}& \rho_{ee}\end{pmatrix}
\text{ which equals }|\psi⟩⟨\psi|
=\begin{pmatrix}C_0\\ C_1\end{pmatrix}\begin{pmatrix}C_0^\ast& C_1^*\end{pmatrix}
\text{ for a pure state}.
$$
Driving the system with Rabi frequency $\Omega$ and detuning $\delta=\omega-\omega_0$, and allowing it to spontaneously emit at a rate $\gamma$ gives the optical Bloch equations:
$$
\begin{align}
\frac{d \rho_{gg}}{dt} & = \gamma \rho_{ee}+\frac i2 (\Omega^* \rho_{eg}-\Omega \rho_{ge}),
\\
\frac{d \rho_{ee}}{dt} & = -\gamma \rho_{ee}+\frac i2 (\Omega \rho_{eg}-\Omega^* \rho_{ge}),
\\
\frac{d \rho_{ge}}{dt} & =-\left(\frac\gamma2+i\delta\right)\rho_{ge}+\frac i2\Omega^*(\rho_{ee}-\rho_{gg}),
\\
\frac{d \rho_{eg}}{dt} & =-\left(\frac\gamma2+i\delta\right)\rho_{eg}+\frac i2\Omega^*(\rho_{gg}-\rho_{ee}).
\end{align}
$$
It is these equations that you should be solving.
