Oscilador armónico cuántico amortiguado: evolución del estado coherente

Estoy tratando de resolver la siguiente ecuación maestra (también similar al oscilador armónico cuántico amortiguado):

$$\frac{d\hat{\rho}}{dt} = \frac{\Gamma}{2}\left(2\hat{a}\hat{\rho}\hat{a}^{\dagger} - \hat{a}^{\dagger}\hat{a} \hat{\rho} - \hat{\rho}\hat{a}^{\dagger}\hat{a} \right)$$ con estado coherente inicial: $\hat{\rho}(0) = |\alpha\rangle\langle\alpha|$. Mi idea era usar la función de Sudarshan y el álgebra D de Gilmore para escribir una ecuación diferencial. El primer paso es descomponer el operador de matriz de densidad en una base de estado coherente: $$\hat{\rho}(t) = \int d^2 \beta\ P(\beta, t) |\beta\rangle \langle \beta|$$ y actuar con los operadores que aparecen en la ecuación inicial:

$$\hat{a}\hat{\rho}\hat{a}^{\dagger} = \int d^2 \beta\ |\beta|^2 P(\beta, t) |\beta\rangle \langle \beta|$$ $$\hat{a}^{\dagger}\hat{a} \hat{\rho} = \int d^2 \beta\ P(\beta,t)\beta\left(\beta^* + \frac{\partial}{\partial\beta} \right)|\beta\rangle \langle \beta| = \int d^2 \beta\ |\beta\rangle\langle\beta| \left(\beta^* - \frac{\partial}{\partial\beta} \right) \beta P(\beta,t)$$ $$\hat{\rho}\hat{a}^{\dagger}\hat{a} = \int d^2 \beta\ P(\beta,t) \beta^*\left(\beta + \frac{\partial}{\partial\beta^*} \right) |\beta\rangle \langle \beta| = \int d^2 \beta\ |\beta\rangle\langle\beta| \left(\beta - \frac{\partial}{\partial\beta^*} \right)\beta^*P(\beta,t)$$ Finalmente obtenemos: $$\frac{\partial P(\beta, t)}{\partial t} = \frac{\Gamma}{2}\left(\beta \frac{\partial}{\partial \beta} + \beta^* \frac{\partial}{\partial \beta^* } +2\right) P(\beta, t)$$

Estoy bastante seguro de que la derivación es correcta porque la ecuación diferencial conserva la traza de la unidad, es decir

$$\frac{d}{dt}\text{Tr}\{\hat{\rho}(t)\} = \int d^2\beta\ \frac{\partial P(\beta,t)}{\partial t} = \frac{\Gamma}{2}\int d^2\beta\ \left(\frac{\partial}{\partial \beta}\beta P(\beta,t) + \frac{\partial}{\partial \beta^*}\beta^* P(\beta,t) \right) = 0$$ Mi idea ahora era usar la exponenciación y escribir la solución casi final como: $$P(\beta, t) = \exp\left[t\frac{\Gamma}{2}\left(\beta \frac{\partial}{\partial \beta} + \beta^* \frac{\partial}{\partial \beta^*} +2\right)\right] \delta^{(2)}(\alpha - \beta)$$ y luego de la definición de la función delta $$\delta^{(2)}(\alpha - \beta) = \frac{1}{\pi^2}\int d^2\eta\ e^{-i\eta^*(\alpha^* - \beta^*)}e^{-i\eta(\alpha - \beta)}$$ puedo escribir $$P(\beta,t) = e^{\Gamma t} \frac{1}{\pi^2} \int d^2\eta\ e^{-i\eta^*\alpha^2} e^{t\frac{\Gamma}{2} \beta^* \partial_{\beta^*}} e^{i\eta^* \beta^*} \times e^{-i\eta \alpha} e^{t\frac{\Gamma}{2} \beta \partial_{\beta}} e^{i\eta \beta} = e^{\Gamma t}\delta^{2}(\alpha - \beta e^{t\Gamma/2})$$ El último paso se puede encontrar aquí . al final me sale $$\rho(t) = e^{\Gamma t}|\alpha e^{-t\Gamma/2}\rangle\langle \alpha e^{-t\Gamma/2}|$$ El problema es que no conserva el rastro a medida que crece hasta el infinito. ¿Alguna idea de dónde cometí un error?