Result

Once we finally take account of all terms, we find that, with $\varepsilon=\pi g(\omega)^2D(\omega)$, and

$$\Delta \omega = P\int^{\infty}_0\frac{g(\omega_j)^2D(\omega_j)}{\omega-\omega_j}d\omega_j$$

$$\frac{d\rho_s}{dt} = -i\Delta \omega[a^{\dagger}a,\rho_s(t)] - \varepsilon(1+\bar{n})\left(\rho_s(t)a^{\dagger}a+a^{\dagger}a\rho_s(t)-2a\rho_s(t)a^{\dagger}\right)$$

$$-\varepsilon(\bar{n})\left(\rho_s(t)aa^{\dagger} + aa^{\dagger}\rho_s(t)-2a^{\dagger}\rho_s(t)a\right)$$