Decay into continuum with complex energies The excited-state amplitude of a single mode coupled to a continuum is governed by the coupled differential equations
$$ \dot c(t) = i\int dk\, g_k^* b_k(t)$$
and 
$$ \dot b_k(t) = -i (\omega_k-\omega_0)b_k(t)+ig_k c(t). $$
Eliminating the continuum degrees of freedom, we arrive at an integro-differential equation for the excited-state amplitude
$$ \dot c(t) =-\left[\int_0^tdt'\int dk\,e^{-i(\omega_k-\omega_0)(t-t')}|g_k|^2\right]c(t'). $$
Often it's permissible to apply the Markov condition $c(t')\to c(t)$. 
In order to perform the resulting integral, a common approach is to approximate the kernel as (at late times)
$$\int_0^t dt'\,e^{-i(\omega_k-\omega_0)(t-t')}\simeq\pi\delta(\omega_k-\omega_0)-i\mathcal P\left(\frac{1}{\omega_k-\omega_0}\right). $$
This only really works when the energies $\omega_k,\omega_0$ are real. Is there a similar simple workaround for complex energies? This arises, e.g., when calculating the decay of a decaying level into another set of modes. A similar problem is considered in https://journals.aps.org/pra/abstract/10.1103/PhysRevA.93.062129, but the formalism seems like overkill for this problem.
 A: Some notes: 


*

*The central equation in the question $$ \dot c(t) =-\left[\int_0^tdt'\int dk\,e^{-i(\omega_k-\omega_0)(t-t')}|g_k|^2\right]c(t') \,, $$ can be solved using Laplace transformations (Krimer2014) without employing the Markov approximation. The method is based on identifying the poles of a certain kernel in the complex spectral plane. For complex $\omega_0$, this would simply shift the poles and should therefore be easy to generalize. Note that the derivation there also uses the Sokhotski–Plemelj theorem alluded to in the question.

*The coupled starting equations are unphysical. Specifically, the interaction picture operators $b_k(t)$ for the field do not fulfill canonical commutation relations and are not associated with physical field degrees of freedom. The reason is the transformation to the interaction picture which is required to be unitary. For real $\omega_0$, this is the case while for complex $\omega_0$, it is not. That said, one can treat a second explicitly Markovian bath that the single mode can decay into by the following equations
$$ \dot c(t) = -i(\omega_0 - i\gamma) c(t) + i\int dk\, g_k^* b_k(t) \,,$$
$$ \dot b_k(t) = -i \omega_kb_k(t)+ig_k c(t)\,, $$
where $\gamma$ is the decay constant/imaginary frequency part. This gives a clear picture of which physical problem we are dealing with and we can get a similar equation to the one in the question. However, this also shows that the manipulation in the question is not particularly convenient. We could for example not transform to an interaction picture, which then gives something like $$ \dot c(t) =-i(\omega_0 - i\gamma) c(t) -\left[\int_0^tdt'\int dk\,e^{-i\omega_k(t-t')}|g_k|^2\right]c(t') \,. $$
Or we couuld transform to a physical interaction picture via $c(t) \rightarrow c(t)e^{-i\omega_0 t}$ to get something like $$ \dot c(t) =-\gamma c(t) -\left[\int_0^tdt'\int dk\,e^{-i(\omega_k-\omega_0)(t-t')}|g_k|^2\right]c(t') \,, $$
where $\omega_0$ is still real. This equation is much easier to interpret.

*In particular, within the Markov approximation ($c(t')\approx c(t)$ over the relevant integration range), this is essentially Wigner-Weisskopf theory of spontaneous emission with an additional atomic decay term. We can now simplify the integrals in exactly the same way as what is done in normal Wigner-Weisskopf theory to get the simple equation $$ \dot c(t) =-(\gamma + \gamma_b) c(t) \,, $$
with the bath decay constant given by $$\gamma_b = \pi |g_\tilde{k}|^2\,.$$ To obtain this, I have used the Sokhotski–Plemelj theorem from the question and $\tilde{k}$ is defined such that $\omega_\tilde{k} = \omega_0$. The $\delta(\omega_k - \omega_0)$-term gives the decay constant. The principal value $\mathcal{P}$-term is technically infinite, but taken to renormalize the bare transition frequency.
As a side note, here is a cool paper (Malekakhlagh2017) that shows how the divergences encountered in the Wigner-Weisskopf theory are unnecessary and a result of starting with equations that are too simple because they do not respect gauge-invariance.
NOTE: In my notation $\gamma$ is minus the imaginary part of $\omega_0$ in the question. My $\omega_0$ is real. So we have $\omega_0^{[\textrm{question}]} = \omega_0 - i\gamma$ in my answer.
