Quantum transport in driven systems with vibrations: Floquet nonequilibrium Green’s functions and the GD approximation

The system is modelled with the following Hamiltonian,

$$H(t)=H_{el}(t)+H_{vib}+H_{e-v}.$$

(1)

The electronic component are given by:

$$\begin{split}H_{el}(t)=H_{central}+H_{leads}(t)+H_{coupling},\end{split}$$

(2)

$$H_{central}(t)=\sum_{ij}\epsilon_{ij}d^{\dagger}_{i}d_{j},$$

(3)

$$H_{leads}(t)=\sum_{k,\alpha=L,R}\epsilon_{k\alpha}(t)c^{\dagger}_{k\alpha}c_{k\alpha},$$

(4)

$$H_{coupling}(t)=\sum_{ik\alpha=L,R}t_{k\alpha i}\;c^{\dagger}_{k\alpha}d_{i}+\mbox{h.c.}$$

(5)

and the Hamiltonian for the vibrations is given as

$$H_{vibratons}=\sum_{\alpha}\omega_{\alpha}a^{\dagger}_{\alpha}a_{\alpha},$$

(6)

and the coupling between the electronic and vibrational components is given by

$$H_{e-v}=\sum_{ij,\alpha}\lambda^{\alpha}_{i,j}Q_{\alpha}d^{\dagger}_{i}d_{j}.$$

(7)

Here, $d_{i}$ ($d^{\dagger}_{i}$) and $c_{k\alpha}$ ($c^{\dagger}_{k\alpha})$ are annihilation (creation) operators for the site $i$ in the central region and $k\alpha$ in the leads, respectively. For the phonons, the quantum operator for position is given by $Q_{\alpha}=\frac{1}{\sqrt{2}}\left(a_{\alpha}+a^{\dagger}_{\alpha}\right)$ and the momentum given by $P_{\alpha}=\frac{1}{\sqrt{2}i}\left(a_{\alpha}-a^{\dagger}_{\alpha}\right)$, where $a_{\alpha}$ $(a^{\dagger}_{\alpha})$ is the annihilation (creation) operator for the phonon $\alpha$. Within the investigation, only the energies of the leads are considered to be time-dependent. It is a simple extension to consider time-dependence within the central energy levels, $\epsilon_{ij}$ and couplings, $t_{k\alpha,i}$.

To capture the dynamics of the system out of equilibrium, we make use of a nonequilibrium Green’s functions approach Schuler2016 ; Haug2008 ; Stefanucci2013 . On the Keldysh contour, we have the electronic contour Green’s function:

$$G_{ij}(\tau,\tau^{\prime})=-i\left\langle T_{c}\left(d_{i}\left(\tau\right)d^{\dagger}_{j}\left(\tau^{\prime}\right)\right)\right\rangle.$$

(8)

Similarly, we have the phononic Green’s functions,

$$\begin{split}D_{\alpha\beta}(\tau,\tau^{\prime})=-i\left\langle T_{c}\left(\Delta Q_{\alpha}\left(\tau\right)\Delta Q_{\beta}\left(\tau^{\prime}\right)\right)\right\rangle\\ =-i\left[\left\langle T_{c}\left(Q_{\alpha}\left(\tau\right)Q_{\beta}\left(\tau^{\prime}\right)\right)\right\rangle-\left\langle Q_{\alpha}\left(\tau\right)\right\rangle\left\langle Q_{\beta}\left(\tau^{\prime}\right)\right\rangle\right]\end{split}$$

(9)

and

$$\begin{split}D^{PP}_{\alpha\beta}(\tau,\tau^{\prime})=-i\left\langle T_{c}\left(\Delta P_{\alpha}\left(\tau\right)\Delta P_{\beta}\left(\tau^{\prime}\right)\right)\right\rangle\\ =-i\left[\left\langle T_{c}\left(P_{\alpha}\left(\tau\right)P_{\beta}\left(\tau^{\prime}\right)\right)\right\rangle-\left\langle P_{\alpha}\left(\tau\right)\right\rangle\left\langle P_{\beta}\left(\tau^{\prime}\right)\right\rangle\right],\end{split}$$

(10)

where $\Delta Q_{\alpha}\left(\tau\right)=Q_{\alpha}\left(\tau\right)-\left\langle Q_{\alpha}\left(\tau\right)\right\rangle$ and $\Delta P_{\alpha}\left(\tau\right)=P_{\alpha}\left(\tau\right)-\left\langle P_{\alpha}\left(\tau\right)\right\rangle$. The corresponding noninteracting Green’s functions are denoted with lower case lettering.

The current from the left or right lead is given by,

$$\begin{split}I_{\alpha}(t)=2Re\left\{\int^{\infty}_{-\infty}dt_{1}Tr\left[G^{<}(t,t_{1})\Sigma_{\alpha}^{A}(t_{1},t)\right.\right.\\ \left.\left.+G^{R}(t,t_{1})\Sigma_{\alpha}^{<}(t_{1},t)\right]\right\},\end{split}$$

(11)

and the occupation of the electrons within the central region is given by

$$n^{el}_{i}(t)=-iG_{ii}^{<}(t,t).$$

(12)

For the phonons, we are principally interested in the phonon occupation:

$$\begin{split}n^{ph}_{\alpha}(t)=\left\langle a^{\dagger}_{\alpha}\left(t\right)a_{\alpha}\left(t\right)\right\rangle\\ =\frac{1}{2}\left[\langle{P_{\alpha}(t)^{2}}\rangle+\langle{Q_{\alpha}(t)^{2}}\rangle\right]-\frac{1}{2}\\ =\frac{1}{2}\left[iD^{<}_{\alpha\alpha}\left(t,t\right)+\left\langle Q_{\alpha}(t)\right\rangle^{2}\right.\\ \left.+iD^{<,PP}_{\alpha\alpha}\left(t,t\right)+\left\langle P_{\alpha}(t)\right\rangle^{2}\right]-\frac{1}{2}.\end{split}$$

(13)

To calculate the electronic and phononic nonequilibrium Green’s functions, we use the Kadanoff-Baym equations:

$$\begin{split}i\frac{\partial}{\partial\tau}G_{ij}\left(\tau,\tau^{\prime}\right)-\sum_{k}\epsilon_{ik}G_{kj}\left(\tau,\tau^{\prime}\right)\\ -\sum_{k}\int_{c}d\tau_{1}\Sigma_{ik}\left(\tau,\tau_{1}\right)G_{kj}\left(\tau_{1},\tau^{\prime}\right)=\delta_{ij}\delta_{c}\left(\tau-\tau^{\prime}\right)\end{split}$$

(14)

and

$$\begin{split}\frac{-1}{\omega_{\alpha}}\left(\frac{d^{2}}{d\tau^{2}}+\omega_{\alpha}^{2}\right)D_{\alpha\beta}\left(\tau,\tau^{\prime}\right)=\delta_{c}\left(\tau-\tau^{\prime}\right)\delta_{\alpha\beta}\\ +\sum_{\gamma}\int d\tau_{1}\Pi_{\alpha\gamma}(\tau,\tau_{1})D_{\gamma\beta}(\tau_{1},\tau^{\prime}).\end{split}$$

(15)

Here, the external influences on the central electrons and phonons is captured in $\Sigma_{ij}(\tau,\tau^{\prime})$ and $\Pi_{\alpha\beta}(\tau,\tau^{\prime})$, respectively. The collective influence on the central regions electrons can be separated into that due to the phonons and the leads:

$$\Sigma\left(\tau,\tau^{\prime}\right)=\Sigma_{int}\left(\tau,\tau^{\prime}\right)+\Sigma_{leads}\left(\tau,\tau^{\prime}\right)$$

(16)

and, similarly for the phonons,

$$\Pi\left(\tau,\tau^{\prime}\right)=\Pi_{int}\left(\tau,\tau^{\prime}\right)+\Pi_{bath}\left(\tau,\tau^{\prime}\right).$$

(17)

Calculating the electronic lead self-energies follows the standard procedure:

$$\Sigma^{<,>,R,A}_{\alpha,ij}(t,t^{\prime})=\sum_{k,k^{\prime}}t^{*}_{k\alpha,i}(t)\;g^{<,>,R,A}_{k\alpha,k^{\prime}\alpha}\left(t,t^{\prime}\right)t_{k^{\prime}\alpha,j}(t^{\prime}),$$

(18)

where the noninteracting lead self-energies follow the standard definitions,

$$g_{k\alpha,k^{\prime}\alpha^{\prime}}^{<}(t,t^{\prime})=if_{k\alpha}e^{-i\int_{t^{\prime}}^{t}dt_{1}\epsilon_{k\alpha}(t_{1})}\delta_{k,k^{\prime}},$$

(19)

$$g_{k\alpha,k^{\prime}\alpha}^{>}(t,t^{\prime})=-i(1-f_{k\alpha})e^{-i\int_{t^{\prime}}^{t}dt_{1}\epsilon_{k\alpha}(t_{1})}\delta_{k,k^{\prime}},$$

(20)

$$g^{R}_{ka,k^{\prime}a}(t,t^{\prime})=-i\Theta\left(t-t^{\prime}\right)e^{-i\int_{t^{\prime}}^{t}dt_{1}\epsilon_{k\alpha}(t_{1})}\delta_{k,k^{\prime}},$$

(21)

$$g^{A}_{ka,k^{\prime}a}(t,t^{\prime})=i\Theta\left(t^{\prime}-t\right)e^{-i\int_{t^{\prime}}^{t}dt_{1}\epsilon_{k\alpha}(t_{1})}\delta_{k,k^{\prime}}.$$

(22)

The time-dependence takes the form $\epsilon_{k\alpha}(t)=\epsilon_{k\alpha}+\phi_{\alpha}(t)$, which allows us to separate out the phase induced by the varying energies of the leads from rest of the self-energy:

$$\begin{split}\Sigma_{\alpha,ij}(t,t^{\prime})=\Sigma^{\prime}_{\alpha,ij}(t-t^{\prime})e^{-i\int_{t^{\prime}}^{t}dt_{1}\phi_{\alpha}(t_{1})}\\ =e^{-i\Phi_{\alpha}(t)}\Sigma^{\prime}_{\alpha,ij}(t-t^{\prime})e^{i\Phi_{\alpha}(t^{\prime})},\end{split}$$

(23)

where $\Phi_{\alpha}(t)$ is the anti-derivative of $\phi_{\alpha}(t)$ and $\Sigma^{\prime}_{\alpha,ij}(t-t^{\prime})$ is the self-energy of the equivalent static case, which are taken in the wide-band approximation:

$$\Sigma^{A/R}_{\alpha,ij}(\omega)=\pm\frac{i}{2}\Gamma_{\alpha,ij},$$

(24)

$$\Sigma^{<}_{\alpha,ij}(\omega)=if_{\alpha}\left(\omega\right)\Gamma_{\alpha,ij},$$

(25)

and

$$\Sigma^{>}_{\alpha,ij}(\omega)=-i\left(1-f_{\alpha}\left(\omega\right)\right)\Gamma_{\alpha,ij},$$

(26)

where the Fermi-Dirac occupation is given by:

$$f_{k\alpha}=\frac{1}{1+e^{(\epsilon_{k\alpha}-\mu_{\alpha})/T_{\alpha}}}.$$

(27)

Within the investigation, the lead energies were driven sinusoidally,

$$\phi_{\alpha}(t)=\Delta_{\alpha}\cos(\Omega_{\alpha}t),$$

(28)

which gives $\Phi(t)=\left(\Delta_{\alpha}/\Omega_{\alpha}\right)\sin\left(\Omega_{\alpha}t\right)$, which can be expressed as a Fourier series with the Jacobi-Anger expansion:

$$e^{i\frac{\Delta_{\alpha}}{\Omega_{\alpha}}\sin(\Omega_{\alpha}t)}=\sum_{n=-\infty}^{n=\infty}J_{n}\left(\frac{\Delta_{\alpha}}{\Omega_{\alpha}}\right)e^{in\Omega_{\alpha}t},$$

(29)

where $J_{n}(x)$ are Bessel functions of the first kind.

For the noninteracting phonons, we have the following phonon Green’s functions:

$$d_{\alpha\beta}^{R}(\omega)=\left[\frac{\frac{1}{2}}{\omega-\omega_{\alpha}+i\eta_{\alpha}}-\frac{\frac{1}{2}}{\omega+\omega_{\alpha}+i\eta_{\alpha}}\right]\delta_{\alpha\beta},$$

(30)

$$d_{\alpha\beta}^{A}(\omega)=\left[\frac{\frac{1}{2}}{\omega-\omega_{\alpha}-i\eta_{\alpha}}-\frac{\frac{1}{2}}{\omega+\omega_{\alpha}-i\eta_{\alpha}}\right]\delta_{\alpha\beta},$$

(31)

$$\begin{split}d_{\alpha\beta}^{<}(\omega)=-\pi i\left[f_{B}(\omega_{\alpha})\delta(\omega-\omega_{\alpha})\right.\\ \left.+(1+f_{B}(\omega_{\alpha}))\delta(\omega+\omega_{\alpha})\right]\delta_{\alpha\beta},\end{split}$$

(32)

$$\begin{split}d_{\alpha\beta}^{>}(\omega)=-\pi i\left[f_{B}(\omega_{\alpha})\delta(\omega+\omega_{\alpha})\right.\\ \left.+(1+f_{B}(\omega_{\alpha}))\delta(\omega-\omega_{\alpha})\right]\delta_{\alpha\beta}.\end{split}$$

(33)

Instead of letting $\eta_{\alpha}$ be infinitesimals, we take them as finite, so to capture the influence of a phonon bath on central phonons. Making use of fluctuation-dissipation relations, we can introduce $\eta_{\alpha}$ into the lesser and greater phonon Green’s functions:

$$\begin{split}d_{\alpha}^{<}(\omega)=\left(d_{\alpha}^{R}(\omega)-d_{\alpha}^{A}\left(\omega\right)\right)f_{B}(\omega)\\ =\left(-\frac{i\eta_{\alpha}}{\left(\omega-\omega_{\alpha}\right)^{2}+\eta_{\alpha}^{2}}+\frac{i\eta_{\alpha}}{\left(\omega+\omega_{\alpha}\right)^{2}+\eta_{\alpha}^{2}}\right)f_{B}(\omega)\end{split}$$

(34)

and

$$\begin{split}d_{\alpha}^{>}(\omega)=\left(d_{\alpha}^{R}(\omega)-d_{\alpha}^{A}\left(\omega\right)\right)\left(1+f_{B}(\omega)\right)\\ =\left(-\frac{i\eta_{\alpha}}{\left(\omega-\omega_{\alpha}\right)^{2}+\eta_{\alpha}^{2}}+\frac{i\eta_{\alpha}}{\left(\omega+\omega_{\alpha}\right)^{2}+\eta_{\alpha}^{2}}\right)\left(1+f_{B}(\omega)\right),\end{split}$$

(35)

where taking the limit of $\eta$,

$$\lim_{\eta\rightarrow 0^{+}}\left[\frac{\eta}{\left(\omega-\omega_{0}\right)^{2}+\eta^{2}}\right]\rightarrow\pi\delta\left(\omega-\omega_{0}\right),$$

(36)

and substituting for $f^{0}_{B}(-\omega)=-\left(f^{0}_{B}(\omega)+1\right)$, gives us back the equations (32) and (33).

For the lesser and greater phonon Green’s functions, we can use the fluctuation-dissipation rules to cast the effects due to the infinitesimals as self-energies. Equating equations (32) and (33) with the associated Keldysh equations gives us:

$$\Pi_{\alpha}^{</>}\left(\omega\right)=\mp 4i\eta_{\alpha}\left(\frac{\omega}{\omega_{\alpha}}\right)f_{B}(\pm\omega).$$

(37)

This phonon self-energy is used to incorporate the effects of the infinitesimal into self-consistent calculations (see equation (17)).

In addition to the above Green’s functions, we need to calculate $\left\langle Q_{\alpha}(\tau)\right\rangle$, $\left\langle P_{\alpha}(\tau)\right\rangle$ and $D^{PP}_{\alpha,\alpha^{\prime}}\left(\tau,\tau^{\prime}\right)$, allowing for the calculation of the phonon occupation. The equation of motion for the average position is given as Schuler2016 ,

$$\begin{split}\frac{-1}{\omega_{\alpha}}\left(\frac{d^{2}}{d\tau^{2}}+\omega_{\alpha}^{2}\right)\left\langle Q_{\alpha}(\tau)\right\rangle=\sum_{ij}-i\lambda^{\alpha}_{ij}G_{ji}\left(\tau,\tau^{+}\right)\\ +\int d\tau_{1}\Pi^{bath}_{\alpha}(\tau,\tau_{1})\left\langle Q_{\alpha}(\tau_{1})\right\rangle.\end{split}$$

(38)

For the terms with momentum operators, we make use of

$$\frac{d}{d\tau}Q_{\alpha}(\tau)=\omega_{\alpha}P_{\alpha}(\tau),$$

(39)

which gives us the equations of motions

$$\frac{d}{d\tau}\langle{Q_{\alpha}(\tau)}\rangle=\omega_{\alpha}\langle{P_{\alpha}(\tau)}\rangle$$

(40)

and

$$\begin{split}\frac{d}{d\tau d\tau^{\prime}}D_{\alpha\beta}\left(\tau,\tau^{\prime}\right)=\omega_{\alpha}\delta_{\alpha\beta}\delta\left(\tau,\tau^{\prime}\right)\\ +D^{PP}_{\alpha\beta}\left(\tau,\tau^{\prime}\right)\omega_{\alpha}\omega_{\beta}.\end{split}$$

(41)

These objects allows us to calculate occupation:

$$\begin{split}n^{ph}_{\alpha}(\tau)=\left\langle T_{c}\left(a_{\alpha}\left(\tau\right)a_{\alpha}^{\dagger}\left(\tau^{+}\right)\right)\right\rangle\\ =\frac{1}{2}\left[\left\langle T_{c}\left(Q^{2}_{\alpha}\left(\tau\right)\right)\right\rangle+\left\langle T_{c}\left(P^{2}_{\alpha}\left(\tau\right)\right)\right\rangle\right]-\frac{1}{2}\\ =\frac{1}{2}\left[iD_{\alpha\alpha}\left(\tau,\tau^{+}\right)+\left\langle Q_{\alpha}(\tau)\right\rangle^{2}\right.\\ \left.+iD^{PP}_{\alpha\alpha}\left(\tau,\tau^{+}\right)+\left\langle P_{\alpha}(\tau)\right\rangle^{2}\right]-\frac{1}{2}.\end{split}$$

(42)

The interaction between electrons and phonons can be approximated within the GD approximation Sakkinen2015 ; Karlsson2020 :

$$\mathbf{\Sigma}_{ij}^{int}(\tau,\tau^{\prime})=\mathbf{\Sigma}_{ij}^{har}(\tau,\tau^{\prime})+\mathbf{\Sigma}_{ij}^{XC}(\tau,\tau^{\prime}),$$

(43)

where

$$\begin{split}\mathbf{\Sigma}_{ij}^{har}(\tau,\tau^{\prime})\\ =-i\delta\left(\tau,\tau^{\prime}\right)\sum_{\beta}\lambda^{\beta}_{ij}\int_{c}d\tau_{1}d_{\beta}(\tau,\tau_{1})\sum_{ml}\lambda_{ml}^{\beta}G_{lm}(\tau_{1},\tau_{1}^{+}),\end{split}$$

(44)

$$\mathbf{\Sigma}^{XC}_{ij}(\tau,\tau^{\prime})=i\sum_{\mu\nu,ml}D_{\mu\nu}\left(\tau,\tau^{\prime}\right)\lambda^{\mu}_{im}G_{ml}\left(\tau,\tau^{\prime}\right)\lambda^{\nu}_{lj},$$

(45)

and

$$\Pi^{int}_{\alpha\beta}(\tau,\tau^{\prime})=-i\sum_{mlkp}\lambda^{\alpha}_{ml}G_{lk}\left(\tau,\tau^{\prime}\right)\lambda^{\beta}_{kp}G_{pm}\left(\tau^{\prime},\tau\right).$$

(46)

Moving the equations of motion from the contour to real time with the greater, lesser, retarded and advanced projections Haug2008 ; Rammer2007 , we can solve the equations of motion with a Floquet approach Honeychurch2020 ; Brandes1997 , where we assume that the system is time-periodic around the central time $T=\frac{t+t^{\prime}}{2}$, and complete a Fourier transform with respect to the relative time, $\tau=t-t^{\prime}$:

$$A(t,t^{\prime})=A(T,\tau)=\sum^{\infty}_{n=-\infty}A(\tau,n)e^{\Omega inT}$$

(47)

and

$$A\left(\omega,n\right)=\frac{1}{P}\int^{P}_{0}dT\;e^{-i\Omega nT}\int^{\infty}_{-\infty}d\tau e^{i\omega\tau}A(T,\tau).$$

(48)

For solving the convolutions of the form

$$C(t,t^{\prime})=\int dt_{1}A(t,t_{1})B(t_{1},t^{\prime}),$$

(49)

we recast the Fourier coefficients into a Floquet matrix,

$$\bar{A}\left(\omega,m,n\right)=A\left(\omega+\frac{\Omega}{2}\left(m+n\right),n-m\right),$$

(50)

which allows us to express the convolution as a matrix multiplication,

$$\bar{C}\left(\omega,m,n\right)=\sum_{r=-\infty}^{\infty}\bar{A}\left(\omega,m,r\right)\bar{B}\left(\omega,r,n\right),$$

(51)

allowing for the Kadanoff-Baym equations to be written as a matrix equation:

$$\begin{split}\left(\omega+\Omega m\right)\bar{G}_{ij}^{R/A}\left(\omega,m,n\right)-\sum_{k}\epsilon_{ik}\bar{G}_{kj}^{R/A}\left(\omega,m,n\right)=\\ \delta_{ij}\delta_{mn}+\sum_{k,r}\bar{\Sigma}_{ik}^{R/A}\left(\omega,m,r\right)\bar{G}_{kj}^{R/A}\left(\omega,r,n\right),\end{split}$$

(52)

$$\begin{split}\bar{G}_{ij}^{<}(\omega,m,n)=\\ \sum_{kw,rs}\bar{G}_{ik}^{R}(\omega,m,r)\bar{\Sigma}^{<}_{kw}\left(\omega,r,s\right)\bar{G}_{wj}^{A}\left(\omega,s,n\right),\end{split}$$

(53)

$$\begin{split}\frac{-1}{\omega_{\alpha}}\left(\omega^{2}_{\alpha}-\left(\omega+\Omega m\right)^{2}\right)\bar{D}^{R}_{\alpha\beta}\left(\omega,m,n\right)=\delta_{\alpha\beta}\delta_{mn}\\ +\sum_{\gamma,r}\bar{\Pi}_{\alpha\gamma}\left(\omega,m,r\right)\bar{D}_{\gamma\beta}\left(\omega,r,n\right),\end{split}$$

(54)

$$\begin{split}\bar{D}_{\alpha\beta}^{<}(\omega,m,n)\\ =\sum_{\nu\gamma,rs}\bar{D}_{\alpha\nu}^{R}(\omega,m,r)\bar{\Pi}^{<}_{\nu\gamma}\left(\omega,r,s\right)\bar{D}_{\gamma\beta}^{A}\left(\omega,s,n\right).\end{split}$$

(55)

Other important objects transform in a similar manner. The lead self-energies, equation (23), making use of equations (29) and (51), transform to

$$\bar{\Sigma}_{\alpha,ij}\left(\omega,m,n\right)=\sum_{pq,lk}\bar{S}_{\alpha}\left(m,p\right){\bar{\Sigma}^{\prime}}_{\alpha,lk}\left(\omega,p,q\right){\bar{S}}_{\alpha}\left(n,q\right),$$

(56)

where $\bar{S}_{\alpha}\left(m,n\right)=J_{m-n}\left(\Delta_{\alpha}/\Omega_{\alpha}\right)$. The Fourier coefficients of the current and occupation can be taken from the Floquet matrices derived from equations (11) and (12):

$$\begin{split}I_{\alpha}(n-m)=\bar{I}_{\alpha}\left(m,n\right)\\ =2\int^{\infty}_{-\infty}\frac{d\omega}{2\pi}\sum_{r,ij}\left[\bar{G}_{ij}^{R}(\omega,m,r)\bar{\Sigma}^{<}_{ji}\left(\omega,r,n\right)\right.\\ \left.+\bar{G}_{ij}^{<}(\omega,m,r)\bar{\Sigma}^{R}_{ji}\left(\omega,r,n\right)\right]\end{split}$$

(57)

and

$$\begin{split}n_{j}^{el}\left(n-m\right)=\bar{n}_{j}^{el}\left(m,n\right)=-i\int^{\infty}_{-\infty}\frac{d\omega}{2\pi}\bar{G}_{jj}^{<}(\omega,m,n),\end{split}$$

(58)

The time-resolved observables can then be found with equation (47), where $t\rightarrow t^{\prime}$, and by taking only the real part, as per equations (12) and (11).

Unfortunately, the interaction self-energies cannot be brought to such an amenable form, and must be calculated from the Fourier coefficients. The interaction self-energies follow the forms

$$C_{XC}(t,t^{\prime})=A\left(t,t^{\prime}\right)B\left(t,t^{\prime}\right),$$

(59)

$$C_{POL}(t,t^{\prime})=A\left(t,t^{\prime}\right)B\left(t^{\prime},t\right),$$

(60)

$$C_{HAR}(t,t^{\prime})=\delta\left(t-t^{\prime}\right)\int^{\infty}_{-\infty}dt_{1}A(t,t_{1})B(t_{1},t_{1}).$$

(61)

Applying the Floquet tranformation to the above gives

$$\begin{split}C_{XC}\left(\omega,n\right)\\ =\sum_{m=-\infty}^{\infty}\int\frac{d\omega^{\prime}}{2\pi}A(\omega^{\prime},m)B(\omega-\omega^{\prime},n-m),\end{split}$$

(62)

$$\begin{split}C_{POL}\left(\omega,n\right)\\ =\sum_{m=-\infty}^{\infty}\int\frac{d\omega^{\prime}}{2\pi}A(\omega^{\prime},m)B(\omega^{\prime}-\omega,n-m)\end{split}$$

(63)

and

$$\begin{split}C_{HAR}(\omega,n)\\ =\sum^{\infty}_{m=-\infty}A\left(-\frac{\Omega}{2}\left(n+m\right),n-m\right)B^{\prime}\left(n\right).\end{split}$$

(64)

The $C_{HAR}\left(\omega,n\right)$ simplifies when $A(t,t^{\prime})=A(t-t^{\prime})$, giving us

$$\begin{split}C_{HAR}\left(\omega,n\right)\\ =\sum^{\infty}_{m=-\infty}\delta_{m,n}A(-\frac{\Omega}{2}\left(m+n\right),0)B(n)\\ =A(-\Omega n,0)B(n).\end{split}$$

(65)

The above allow us to cast the interaction self-energies in terms of their Fourier coefficients:

$$\begin{split}\Sigma^{R}_{har,ij}\left(\omega,r\right)=\\ -i\sum_{\beta}\lambda^{\beta}_{ij}d^{R}_{\beta}(\omega=-\Omega r)\int\frac{d\omega}{2\pi}\sum_{kw}\lambda_{kw}^{\beta}G_{wk}^{<}(\omega,r)\end{split}$$

(66)

$$\begin{split}\Sigma^{R}_{XC,ij}\left(\omega,r\right)=\int\frac{d\omega^{\prime}}{2\pi}\sum^{\infty}_{n=-\infty}\sum_{\mu\nu,ml}\\ iD^{<}_{\mu\nu}\left(\omega^{\prime},n\right)\lambda^{\mu}_{im}G^{R}_{ml}\left(\omega-\omega^{\prime},r-n\right)\lambda^{\nu}_{lj}\\ +iD^{R}_{\mu\nu}\left(\omega^{\prime},n\right)\lambda^{\mu}_{im}G^{<}_{ml}\left(\omega-\omega^{\prime},r-n\right)\lambda^{\nu}_{lj}\\ +iD^{R}_{\mu\nu}\left(\omega^{\prime},n\right)\lambda^{\mu}_{im}G^{R}_{ml}\left(\omega-\omega^{\prime},r-n\right)\lambda^{\nu}_{lj}\end{split}$$

(67)

$$\begin{split}\Sigma^{<}_{XC,ij}\left(\omega,r\right)=\int\frac{d\omega^{\prime}}{2\pi}\sum^{\infty}_{n=-\infty}\sum_{\mu\nu,ml}\\ iD^{<}_{\mu\nu}\left(\omega^{\prime},n\right)\lambda^{\mu}_{im}G^{<}_{ml}\left(\omega-\omega^{\prime},r-n\right)\lambda^{\nu}_{lj}\end{split}$$

(68)

$$\begin{split}\Sigma^{>}_{XC,ij}\left(\omega,r\right)=\int\frac{d\omega^{\prime}}{2\pi}\sum^{\infty}_{n=-\infty}\sum_{\mu\nu,ml}\\ iD^{>}_{\mu\nu}\left(\omega^{\prime},n\right)\lambda^{\mu}_{im}G^{>}_{ml}\left(\omega-\omega^{\prime},r-n\right)\lambda^{\nu}_{lj}\end{split}$$

(69)

$$\begin{split}\Pi^{R}_{\alpha\beta}(\omega,r)=\int\frac{d\omega^{\prime}}{2\pi}\sum^{\infty}_{n=-\infty}\sum_{mlkp}\\ -i\lambda^{\alpha}_{ml}G^{<}_{lk}\left(\omega^{\prime},n\right)\lambda^{\beta}_{kp}G^{A}_{pm}\left(\omega^{\prime}-\omega,r-n\right)\\ -i\lambda^{\alpha}_{ml}G^{R}_{lk}\left(\omega^{\prime},n\right)\lambda^{\beta}_{kp}G^{<}_{pm}\left(\omega^{\prime}-\omega,r-n\right)\end{split}$$

(70)

$$\begin{split}\Pi_{\alpha\beta}^{<}(\omega,r)=\int\frac{d\omega^{\prime}}{2\pi}\sum^{\infty}_{n=-\infty}\sum_{mlkp}\\ -i\lambda^{\alpha}_{ml}G^{<}_{lk}\left(\omega^{\prime},n\right)\lambda^{\beta}_{kp}G^{>}_{pm}\left(\omega^{\prime}-\omega,r-n\right)\end{split}$$

(71)

$$\begin{split}\Pi_{\alpha\beta}^{>}(\omega,r)=\int\frac{d\omega^{\prime}}{2\pi}\sum^{\infty}_{n=-\infty}\sum_{mlkp}\\ -i\lambda^{\alpha}_{ml}G^{>}_{lk}\left(\omega^{\prime},n\right)\lambda^{\beta}_{kp}G^{<}_{pm}\left(\omega^{\prime}-\omega,r-n\right).\end{split}$$

(72)

Solving equation (38), we can cast the average phonon positions, and the subsequent average phonon momenta, in terms of Fourier coefficients

$$\begin{split}\left\langle Q_{\alpha}\right\rangle(r)\\ =-id^{R}_{\alpha}\left(\omega=-\Omega r\right)\int\frac{d\omega}{2\pi}\sum_{ml}\lambda_{ml}^{\alpha}G_{lm}^{<}\left(\omega,r\right),\end{split}$$

(73)

$$\left\langle P_{\alpha}\right\rangle(r)=\frac{ir\Omega}{\omega_{0}}\left\langle Q_{\alpha}\right\rangle(r).$$

(74)

We can complete a similar process for the phonon momentum Green’s functions, allowing us to calculate variance of the momentum operators:

$$\begin{split}\left\langle\left(\Delta P_{\alpha}\right)^{2}\right\rangle(r)=iD_{\alpha\alpha}^{PP,<}(r)\\ =\frac{i}{{\omega_{0}}^{2}}\int\frac{d\omega}{2\pi}\left(\omega^{2}-\left(\frac{r\Omega}{2}\right)^{2}\right)D_{\alpha\alpha}^{<}(\omega,r).\end{split}$$

(75)