Model: spin-boson

To derive Fermi’s golden rule rate for a nonequilibrium initial state, we use polaron transformation to remove the diagonal linear vibronic couplings (LVC). We start from the usual LVC Hamiltonian
\begin{align} H_{0} & =\left(E_{1}+\omega a^{\dagger}a+c_{1}(a^{\dagger}+a)\right)\left|1\right\rangle \left\langle 1\right|\\ & +\left(E_{2}+\omega a^{\dagger}a+c_{2}(a^{\dagger}+a)\right)\left|2\right\rangle \left\langle 2\right|\\ V & =g_{12}\left|1\right\rangle \left\langle 2\right|+h.c. \end{align}
Here, only a single vibrational mode is used for simplicity of derivation. To include multiple modes, one only needs to replace the single mode term in the conclusion with a sum over all modes.

Polaron transformation

The polaron transformation uses the unitary operator
\begin{align} U_{p}=e^{\sum_{n}(a^{\dagger}-a)\frac{c_{n}}{\omega}\left|n\right\rangle \left\langle n\right|}. \end{align}

Notice $e^{\lambda(a^{\dagger}-a)}ae^{-\lambda(a^{\dagger}-a)}=a-\lambda$ and apply $U_{p}$ to $H_{0}$:
\begin{align} \tilde{H}_{0}=&U_{p}H_{0}U_{p}^{\dagger}\\ =& \left(E_{1}+\omega(a^{\dagger}-\frac{c_{1}}{\omega})(a-\frac{c_{1}}{\omega})+c_{1}(a^{\dagger}+a-2\frac{c_{1}}{\omega})\right)\left|1\right\rangle \left\langle 1\right|\\ +& \left(E_{2}+\omega a^{\dagger}a-\frac{c_{3}^{2}}{\omega}\right)\left|2\right\rangle \left\langle 2\right|\\ =& \sum_{n}\left(E_{n}+\omega a^{\dagger}a-c_{n}^{2}/\omega\right)\left|n\right\rangle \left\langle n\right|. \end{align}

Apply $U_{p}$ to $V$:
\begin{align} \tilde{V}= & U_{p}VU_{p}^{\dagger}\\ = & g_{12}\left|1\right\rangle \left\langle 2\right|e^{(a^{\dagger}-a)(c_{1}-c_{2})/\omega}+h.c. \end{align}
To use the perturbation theory, first note
\begin{align} \tilde{V}_{I}(t) & =e^{i\tilde{H}_{0}t}\tilde{V}e^{-i\tilde{H}_{0}t}\\ & =g_{12}\left|1\right\rangle \left\langle 2\right|e^{i\left((E_{1}-\frac{c_{1}^{2}}{\omega})-(E_{2}-\frac{c_{2}^{2}}{\omega})\right)t}e^{(a^{\dagger}e^{i\omega t}-ae^{-i\omega t})(c_{1}-c_{2})/\omega}\\ & +g_{12}^{*}\left|2\right\rangle \left\langle 1\right|e^{-i\left((E_{1}-\frac{c_{1}^{2}}{\omega})-(E_{2}-\frac{c_{2}^{2}}{\omega})\right)t}e^{-(a^{\dagger}e^{i\omega t}-ae^{-i\omega t})(c_{1}-c_{2})/\omega} \end{align}

Polulation of state $\left|1\right\rangle $

Then we consider, the population of state $\left|1\right\rangle $ and let the initial nuclear density matrix be $\rho_{n}$:
\begin{align} P_{1}(t) & =\mathrm{Tr}\left[\rho_{n}\left|1\right\rangle \left\langle 1\right|\underbrace{\underbrace{U_{p}^{\dagger}\hat{U}_{I}^{\dagger}(t)e^{iH_{0}t}U_{p}}_{U_{S}^{\dagger}(t)}\left|1\right\rangle \left\langle 1\right|\underbrace{U_{p}^{\dagger}e^{-iH_{0}t}\hat{U}_{I}(t)U_{p}}_{U_{S}(t)}}_{\left|1\right\rangle \left\langle 1\right|\text{in the Herisenberg picture}}\right]\\ & =\mathrm{Tr}_{n}\left[\rho_{n}\left(e^{-(a^{\dagger}-a)\frac{c_{1}}{\omega}}\left\langle 1\right|\hat{U}_{I}^{\dagger}(t)\left|1\right\rangle e^{(a^{\dagger}-a)\frac{c_{1}}{\omega}}\right)\left(e^{-(a^{\dagger}-a)\frac{c_{1}}{\omega}}\left\langle 1\right|\hat{U}_{I}^{\dagger}(t)\left|1\right\rangle e^{(a^{\dagger}-a)\frac{c_{1}}{\omega}}\right)\right]\\ \\ & \approx1-2\mathrm{Re}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\mathrm{Tr}_{n}\left[\rho_{n}e^{-(a^{\dagger}-a)\frac{c_{1}}{\omega}}\left\langle 1\right|\tilde{V}_{I}(t_{1})\tilde{V}_{I}(t_{2})\left|1\right\rangle e^{(a^{\dagger}-a)\frac{c_{1}}{\omega}}\right] \end{align}

Calculation of the integrand

We consider the quantity
\begin{align*} & \left\langle 1\right|\tilde{V}(t_{1})\tilde{V}(t_{2})\left|1\right\rangle \\ = & g_{12}^{2}e^{i\Delta E(t_{1}-t_{2})}\exp\left[(a^{\dagger}e^{i\omega t_{1}}-ae^{-i\omega t_{1}})(c_{1}-c_{2})/\omega\right]\exp\left[-(a^{\dagger}e^{i\omega t_{2}}-ae^{-i\omega t_{2}})(c_{1}-c_{2})/\omega\right] \end{align*}
where $\Delta E=(E_{1}-\frac{c_{1}^{2}}{\omega})-(E_{2}-\frac{c_{2}^{2}}{\omega})$.

Before averaging over the initial nuclear condition $\rho_{n}=e^{-\beta\omega(a^{\dagger}+c_{0})(a+c_{0})}/\mathcal{Z}=e^{-(a^{\dagger}-a)\frac{c_{0}}{\omega}}\frac{e^{-\beta\omega a^{\dagger}a}}{\mathcal{Z}}e^{(a^{\dagger}-a)\frac{c_{0}}{\omega}}$ we need to calculate
\begin{align} & \exp\left[(a^{\dagger}e^{i\omega t_{1}}-ae^{-i\omega t_{1}})(c_{1}-c_{2})/\omega\right]\exp\left[-(a^{\dagger}e^{i\omega t_{2}}-ae^{-i\omega t_{2}})(c_{1}-c_{2})/\omega\right]\\ = & \exp\left[\frac{c_{1}-c_{2}}{\omega}\left(a^{\dagger}(e^{i\omega t_{1}}-e^{i\omega t_{2}})-a(e^{-i\omega t_{1}}-e^{-i\omega t_{2}})\right)\right]\\ \times & \exp\left(\frac{(c_{1}-c_{2})^{2}}{\omega^{2}}\left[a^{\dagger}e^{i\omega t_{1}}-ae^{-i\omega t_{1}},-(a^{\dagger}e^{i\omega t_{2}}-ae^{-i\omega t_{2}})\right]/2\right)\\ = & \exp\left[\frac{c_{1}-c_{2}}{\omega}\left(a^{\dagger}(e^{i\omega t_{1}}-e^{i\omega t_{2}})-a(e^{-i\omega t_{1}}-e^{-i\omega t_{2}})\right)\right]e^{-i\frac{(c_{1}-c_{2})^{2}}{\omega^{2}}\sin\omega(t_{1}-t_{2})} \end{align}
and rotate it by sandwiching it with
\begin{align} & e^{(a^{\dagger}-a)\frac{c_{0}}{\omega}}e^{-(a^{\dagger}-a)\frac{c_{1}}{\omega}}(\cdot)e^{(a^{\dagger}-a)\frac{c_{1}}{\omega}}e^{-(a^{\dagger}-a)\frac{c_{0}}{\omega}}\\ = & e^{-(a^{\dagger}-a)\frac{c_{1}-c_{0}}{\omega}}(\cdot)e^{(a^{\dagger}-a)\frac{c_{1}-c_{0}}{\omega}}. \end{align}

After rotation, the operator expression above is
\begin{align} & e^{-(a^{\dagger}-a)\frac{c_{1}-c_{0}}{\omega}}\exp\left[\frac{c_{1}-c_{2}}{\omega}\left(a^{\dagger}(e^{i\omega t_{1}}-e^{i\omega t_{2}})-a(e^{-i\omega t_{1}}-e^{-i\omega t_{2}})\right)\right]e^{(a^{\dagger}-a)\frac{c_{1}-c_{0}}{\omega}}\\ = & \exp\left[\frac{c_{1}-c_{2}}{\omega}\left((a^{\dagger}+\frac{c_{1}-c_{0}}{\omega})(e^{i\omega t_{1}}-e^{i\omega t_{2}})-(a+\frac{c_{1}-c_{0}}{\omega})(e^{-i\omega t_{1}}-e^{-i\omega t_{2}})\right)\right]. \end{align}
The average value is
\begin{align} & \mathrm{Tr}_{n}\left[\rho_{n}e^{-(a^{\dagger}-a)\frac{c_{1}-c_{0}}{\omega}}\left\langle 1\right|\tilde{V}_{I}(t_{1})\tilde{V}_{I}(t_{2})\left|1\right\rangle e^{(a^{\dagger}-a)\frac{c_{1}-c_{0}}{\omega}}\right]/e^{i\frac{(c_{1}-c_{2})^{2}}{\omega^{2}}\sin\omega(t_{1}-t_{2})}\\ = & \mathrm{Tr}_{n}\left\{ \rho_{n}\exp\left[\frac{c_{1}-c_{2}}{\omega}\left((a^{\dagger}+\frac{c_{1}-c_{0}}{\omega})(e^{i\omega t_{1}}-e^{i\omega t_{2}})-(a+\frac{c_{1}-c_{0}}{\omega})(e^{-i\omega t_{1}}-e^{-i\omega t_{2}})\right)\right]\right\} /e^{i\frac{(c_{1}-c_{2})^{2}}{\omega^{2}}\sin\omega(t_{1}-t_{2})}\\ = & \mathrm{Tr}_{n}\left\{ \rho_{n}\exp\left[\frac{c_{1}-c_{2}}{\omega}\left(a^{\dagger}(e^{i\omega t_{1}}-e^{i\omega t_{2}})-a(e^{-i\omega t_{1}}-e^{-i\omega t_{2}})\right)+\frac{c_{1}-c_{2}}{\omega}\frac{c_{1}-c_{0}}{\omega}\left(2i\sin\omega t_{1}-2i\sin\omega t_{2}\right)\right]\right\} /e^{i\frac{(c_{1}-c_{2})^{2}}{\omega^{2}}\sin\omega(t_{1}-t_{2})}\\ = & e^{\frac{-(c_{1}-c_{2})^{2}}{\omega^{2}}\left[\left(\coth(\beta\omega/2)\right)\left(1-\cos\omega(t_{1}-t_{2})\right)+i\sin\omega(t_{1}-t_{2})\right]}e^{2i\frac{(c_{1}-c_{2})(c_{1}-c_{0})}{\omega^{2}}(\sin\omega t_{1}-\sin\omega t_{2})} \end{align}
where we used $\left\langle e^{\alpha_{1}\hat{a}+\beta_{1}\hat{a}^{\dagger}}e^{\alpha_{2}\hat{a}+\beta_{2}\hat{a}^{\dagger}}\right\rangle _{T}=e^{\left(\alpha_{1}+\alpha_{2}\right)\left(\beta_{1}+\beta_{2}\right)(n+1/2)+(1/2)\left(\alpha_{1}\beta_{2}-\beta_{1}\alpha_{2}\right)}$ and $\langle\cdots\rangle_{T}$ is the thermal average over $\rho_{n}=\frac{e^{-\beta\omega a^{\dagger}a}}{\mathcal{Z}}$ [see Nitzan’s Chemical Dynamics Eq. (12.53)] ¹.

References