text:
\section{Introduction} \label{sec:intro}
Pump-probe spectroscopy is one of the most important tools for studying complex materials with strongly correlated degrees of freedom. The femtosecond time-resolution allows one to examine excitation and relaxation phenomena from different degrees of freedom (electronic, spin, lattice) at their natural timescales. The technique has been used to probe many different materials in the time-domain: the classification of charge-density-wave insulators \cite{CDWinsulators,CDWPumpProbe} and charge-stripe materials \cite{epcOscillations}; examining the charge gap of spin-density-wave materials \cite{SDWPumpProbe}; determining the electron-phonon coupling in a metal \cite{EPCinAmetal}; studying the response in high $T_c$ cuprates \cite{Tc1, Tc2,Tc3,PDHSSC} and in iron-based superconductors \cite{iron}. In a recent experimental study \cite{iron}, evidence of bosonic mode coupling was seen in LiFeAs. The high $T_c$ cuprates are also relevant systems, since the degree of electron-phonon coupling involved in superconductivity in these materials is not completely understood. For example, a recent study \cite{2020eph} indicates the superconducting behavior cannot be understood in terms of electron-electron or electron-phonon coupling alone, but only by considering both effects.

Another class of materials that this theory might be relevant for are those that have been observed to have polarons (or possibly bipolarons) generated by nonequilibrium pumps. It is difficult to experimentally differentiate between polarons and bipolarons, so some of these systems may be creating bipolarons, which naturally occur in the atomic limit and often are strongly localized. In Ref.~\onlinecite{Polaron1}, the authors use pump-probe spectroscopy to measure polaron lifetimes in a mixed-valence perovskite material under high pressures. They also argue that photoexcitation is analogous to applying a high pressure to this system, evidenced by the response of the phonons and the fact that they both lead to an intervalence charge transfer. More recently, experiments show a photoinduced a transition to a long-lived polaronic state in lead halide perovskite \cite{Polaron2}, and in Ref.~\onlinecite{Polaron3} the authors demonstrate a technique for qualitatively estimating the magnitude and shape of a photoinduced polaronic distortion via femtosecond diffuse x-ray scattering techniques.

Improving the capabilities of pump-probe experiments is currently an active field of research \cite{2019techniqueProbe,3ShotPP2020, PPSpec}, and as time and energy resolution increase, we can expect to obtain a clearer picture of the behavior of electron-phonon coupling, electron-electron coupling, and the interplay between them, in complex materials and in nonequilibrium situations.

The pump is an ultrafast and high-power electric field pulse, which may selectively excite the electronic or lattice degrees of freedom in a material. The subsequent probe is often a second electric field pulse, which reveals information about the material after it has been excited by the initial pulse. By varying the time delay between the pump and probe pulses, experiments can ``watch'' as a material relaxes back to equilibrium. A common feature of these experiments is that the transient states exhibit properties not found in equilibrium.

Another interesting question is whether the pump pulse dynamically modifies either the electron-phonon coupling or the electron-electron coupling in pump-probe experiments. The Shen group  demonstrated\cite{EPCEvidence} that the electron-phonon coupling in the high $T_c$ cuprates is related to a kink in the electronic dispersion near the Fermi energy. More recent work concludes that the strength of the kink is determined by a combination of electron-phonon and electron-electron coupling \cite{huang2020KINK}. A softening of the kink has been observed in pump-probe experments\cite{EPCquenching}, but such softening alone is not sufficient to establish a dynamic change in the electron-phonon coupling due to nonequilibrium effects \cite{JimLexKink}. A recent experimental study of electron thermalization in laser-excited nickel\cite{RecentNickelPaper}, suggests a reduction of electron repulsion may explain the observation of an increasing spectral red shift with increasing electron temperature.

The ability to directly extract the electron-phonon coupling from experiment is an area of active interest that has seen exciting recent developments. Namely, a new analysis technique called frequency-domain angle-resolved photo-emission spectroscopy (FD-ARPES) enables observation of the electron-phonon coupling directly from experimental data \cite{DirectEPC}, and has been used to measure the electron-phonon coupling during a photo-induced insulator to metal transition \cite{DirectEPCexperiment}.

In this work, we allow for arbitrary changes in both the electron-electron and electron-phonon coupling without determining their time-dependence from a self-consistent reaction to the pump pulse. Similar studies have been done previously \cite{Werner_2013,Werner_2015} for the lattice model, but we are not aware of work that varies both electron-electron and electron-phonon couplings as functions of time. We also note that the single and two particle atomic Green's functions presented here are the required inputs to a systematic strong-coupling-expansion on the lattice, generalizing previous work \cite{JimSCE}, which applied only to the Hubbard model.
Recently,\cite{SCE2020} a strong-coupling-expansion of the Holstein-Hubbard model with constant coupling parameters found a rich landscape of phases, including distinct pair-density-waves, spin-density-waves, charge-density-waves, and superconducting regions as well as regions of phase separation.

In this work, we provide an exact (equilibrium and nonequilibrium) solution to the atomic Hubbard-Holstein model for arbitrary electron-electron and electron-phonon couplings. In addition to serving as the starting point for strong-coupling-based perturbation theories, it also provides an important benchmark for nonequilibrium calculations, where exact results are rare. The reason why this problem can be solved exactly lies in the fact that the electron number for each spin is conserved in the atomic limit. But, even with that, the dynamics is sufficiently complex that one needs to use advanced techniques to factorize the evolution operator in order to determine the exact results presented here.

The summary of the paper is as follows. In Sec.~\ref{sec:formalism}, we define the time-dependent Hamiltonian and provide general expressions for the spectral moments, which are analyzed in Sec~\ref{sec:moms}. In Sec.~\ref{sec:timedep} we evaluate the time evolution operator, the partition function, and the time dependence of the field operators. These allow us to calculate the single-particle Green's function in Sec.~\ref{sec:onepart}, verify the moments in Sec.~\ref{sec:moms2} and calculate the two-particle Green's function in Sec.~\ref{sec:twopart}. Our conclusions are presented in Sec.~\ref{sec:conclusion}. Appendices summarize some of the more cumbersome formulas.
\section{Formalism} \label{sec:formalism}
 The Holstein-Hubbard model is the simplest model Hamiltonian which describes both electron-electron and electron-phonon interactions \cite{Hubbard,Holstein}. It includes electrons that hop from one lattice site to another and Einstein phonons with the same frequency on each site. The electrons interact with the phonons via a direct coupling of the charge to the phonon coordinate. The electrons also repel each other via an on-site Coulomb repulsion. In the atomic limit, the electrons do not hop. This means we cannot represent electric-field effects with vector potentials and these systems never have electric current flow. The Holstein-Hubbard model for the atomic site then becomes
 \begin{equation}
    \hat{\mathcal H}(t)=\hbar\omega\left (\hat{a}^\dagger\hat{a}+\frac12\right )+\big(g(t)\hat{x}-\mu\big)(\hat{n}_{\uparrow}+\hat{n}_{\downarrow})+U(t)\hat{n}_{\uparrow}\hat{n}_{\downarrow}
\end{equation}
where $\hat{n}_{\sigma}=\hat{c}^{\dagger}_{\sigma}\hat{c}_{\sigma}$ is the fermionic number operator, $\hat{c}_{\sigma}^{\dagger}~(\hat{c}_{\sigma})$ create (destroy) an electron with spin $\sigma$, $\mu$ is the chemical potential, $\hat{a}^{\dagger}~(\hat{a})$ are the bosonic creation (annihilation operators), $\hat{x}=\sqrt{\frac{\hbar}{2m\omega}}(\hat{a}+\hat{a}^{\dagger})$ is the phonon coordinate operator (with $\omega$ the phonon frequency), $g(t)$ is the electron-phonon coupling, and $U(t)$ is the electron-electron coupling. We drop the zero point energy from the Hamiltonian because it just adds an overall constant shift to all energies and we set $\hbar=1$ for the remainder of this work. Note that because $[\hat{\mathcal H}(t),\hat{n}_{\sigma}]=0$, the total electron number (and the electron number for each spin) is conserved. This symmetry is what allows us to derive an exact expression for the retarded Green's function; it is immediately broken once hopping between lattice sites is introduced. Note that we start our system in equilibrium in the limit as $t\to-\infty$ and all of the time dependence in the coupling constants starts at some finite time after $t_{min}$.

The retarded Green's function is defined as
\begin{equation}
\label{gfdef}
    g^R_{\sigma}(t_1,t_2)=-i\theta(t_1-t_2)\big\langle \big\{\hat{c}_{\sigma}(t_1),\hat{c}_{\sigma}^{\dagger}(t_2)\big\}_+\rangle
\end{equation}
where $\theta(t)$ is
the Heaviside function, $\langle \dots\rangle$ indicates the thermal average at the initial time $\frac{1}{\mathcal Z}\text{Tr}\big\{e^{-\beta\hat{\mathcal H}(t_{min})} \dots\big\}$, and all time-dependence of operators is in the Heisenberg representation,
\begin{equation} \hat{O}_H(t)=\hat{U}^{\dagger}(t,t_{min})\hat{O}\hat{U}(t,t_{min})\end{equation} with $\hat{U}(t,t_{min})$ the time-evolution operator, \begin{equation}
    \hat{U}(t,t_{min}) = \mathcal{T}_t\exp\left[-i\int_{t_{min}}^t\!\!\!dt'\ \hat{\mathcal H}(t')\right].
\end{equation}
Here, $\mathcal{T}_t$ is the time-ordering operator, which places the latest times to the left and the Hamiltonian is in the Schr\" odinger representation; $t_{min}$ is a time where the system is in thermal equilibrium and is before any coupling constants start varying with time. The spectral moments of the Green's function are defined via
\begin{equation}
    \mu^{Rn}_{\sigma}(t_{ave})=-\frac{1}{\pi}\int_{-\infty}^{\infty}\,d\omega\ \omega^n\text{Im}\,g^R_{\sigma}(t_{ave},\omega)
\end{equation}
where
\begin{eqnarray}
    &~&g_{\sigma}^R(t_{ave},\omega)=\nonumber\\
    &~&~~\int_0^{\infty}\,dt_{rel}\ e^{i\omega t_{rel}}g^R_{\sigma}\left (t_{ave}+\frac{t_{rel}}{2},t_{ave}-\frac{t_{rel}}{2}\right )
\end{eqnarray}
and we are now using Wigner coordinates of average and relative time,  defined by $t_{ave}=\frac{t_{1}+t_{2}}{2}$ and $t_{rel}=t_1-t_2$. The moments can be rewritten as derivatives with respect to relative time via
\begin{align}
\label{derivs}
    &\mu_{\sigma}^{Rn}(t_{ave})=\lim_{t_{rel}\to 0}\\
    &~~~\text{Re}\bigg\langle i^{n}\frac{d^n}{dt^n_{rel}}\left\{\hat{c}_{\sigma}\left (t_{ave}+\frac{t_{rel}}{2}\right ),\hat{c}_{\sigma}^{\dagger}\left (t_{ave}-\frac{t_{rel}}{2}\right )\right\}_+\bigg\rangle.\nonumber
\end{align}
 Finally we define the time-ordered two-particle atomic Green's function by
 \begin{equation}
 \label{2partgf}
    \mathcal{G}_{\sigma\sigma'}(t_0,t_1,t_2,t_3)=-\langle \mathcal{T}_t \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma'}(t_1)\hat{c}_{\sigma'}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle.
\end{equation}
We will solve for both of these correlation functions below.
\section{Spectral Moment Sum Rules} \label{sec:moms}
 By repeated application of Heisenberg's equation of motion, the spectral moment sum rules can be written as commutators of the field operators with the Hamiltonian plus terms that arise from derivatives of time-dependent coupling constants in the Hamiltonian. After some extensive algebra, we find that the first four moments satisfy \cite{moment2,moments1,erratum,erratum2}
\begin{equation}
    \mu^{R0}_{\sigma}(t_{ave})=\text{Re}\bigg\langle\big\{\hat{c}_{\sigma}(t_{ave}),\hat{c}_{\sigma}^{\dagger}(t_{ave})\big\}_+\bigg\rangle,
\end{equation}
\begin{equation}
    \mu_{\sigma}^{R1}(t_{ave})=\text{Re}\bigg\langle \big\{\big[\hat{c}_{\sigma}(t_{ave}),\hat{\mathcal H}(t_{ave})\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\big\}_+\bigg\rangle,
\end{equation}
\begin{widetext}
\begin{equation}
    \mu_{\sigma}^{R2}(t_{ave})=\text{Re}\bigg\langle \bigg\{\big[[\hat{c}_{\sigma}(t_{ave}),\hat{\mathcal H}(t_{ave})],\hat{\mathcal H}(t_{ave})\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle,
\end{equation}
and
\begin{align}
\begin{split}
    \mu_{\sigma}^{R3}(t_{ave})=\text{Re}\bigg\langle \bigg\{\bigg[\big[[\hat{c}_{\sigma}(t_{ave}),\hat{\mathcal H}(t_{ave})],\hat{\mathcal H}(t_{ave})\big],\hat{\mathcal H}(t_{ave})\bigg],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle\\ +\text{Re}\frac{i}{2}\bigg(\bigg\langle\bigg\{\big[[\hat{c}_{\sigma}(t_{ave}),\frac{\partial \hat{\mathcal H}(t_{ave})}{\partial t_{ave}}],\hat{\mathcal H}(t_{ave})\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle +\bigg\langle\bigg\{\big[[\hat{c}_{\sigma}(t_{ave}),\hat{\mathcal H}(t_{ave})],\frac{\partial \hat{\mathcal H}(t_{ave})}{\partial t_{ave}}\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle\bigg)\\ -\frac{1}{4}\text{Re}\bigg\langle\big[\hat{c}_{\sigma}(t_{ave}),\frac{\partial^2 \hat{\mathcal H}(t_{ave})}{\partial t_{ave}^2}\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle.
    \end{split}
\end{align}
\end{widetext}
Evaluating the commutators in the atomic limit, we obtain the following:
\begin{equation}
    \mu_{\sigma}^{R0}(t_{ave})=1,
\end{equation}
\begin{equation}
\label{mom1}
    \mu_{\sigma}^{R1}=-\big(\mu-U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle -g(t_{ave})\langle \hat{x}(t_{ave})\rangle\big),
\end{equation}
\begin{align}
\label{mom2}
    \mu_{\sigma}^{R2}(t_{ave})&=\big(\mu-U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle -g(t_{ave})\langle \hat{x}(t_{ave})\rangle\big)^2\\
    &+U(t_{ave})^2
    \big[\langle \hat{n}_{\bar{\sigma}}\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^2\big]\nonumber\\
    &+g^2(t_{ave})\big[\langle \hat{x}^2(t_{ave})\rangle-\langle \hat{x}(t_{ave})\rangle^2\big]\nonumber\\
    &+2U(t_{ave})g(t_{ave})\big[\langle \hat{n}_{\bar{\sigma}}\hat{x}(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle \langle \hat{x}(t_{ave})\rangle\big],\nonumber
\end{align}
\begin{widetext}
\begin{align}
\label{mom3}
    \mu_{\sigma}^{R3}(t_{ave})&=-\big(\mu-g(t_{ave})\langle \hat{x}(t_{ave})\rangle -U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle\big)^3+g^3(t_{ave})\big(\langle \hat{x}^3(t_{ave})\rangle -\langle \hat{x}(t_{ave})\rangle^3\big)\\&-3\mu g^2(t_{ave})\big(\langle \hat{x}^2(t_{ave})\rangle-\langle \hat{x}(t_{ave})\rangle^2\big)
    +U^3(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^3\big)-3\mu U^2(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^2\big)\nonumber\\
    &-6\mu g(t_{ave})U(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\hat{x}(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle\langle \hat{x}(t_{ave})\rangle\big)+3U^2(t_{ave})g(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\hat{x}(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^2\langle \hat{x}(t_{ave})\rangle\big)\nonumber\\
    &+3U(t_{ave})g^2(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\hat{x}^2(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle\langle \hat{x}(t_{ave})\rangle^2 \big)\nonumber\\
    &+g(t_{ave})\omega^2\langle \hat{x}(t_{ave})\rangle +\frac{g^2(t_{ave})}{2m}(2\langle \hat{n}_{\bar{\sigma}}\rangle +1) -\frac{1}{4}\text{Re}\bigg(\frac{d^2g(t_{ave})}{dt_{ave}^2}\langle \hat{x}(t_{ave})\rangle +\frac{d^2U(t_{ave})}{dt_{ave}^2}\langle \hat{n}_{\bar{\sigma}}\rangle \bigg)\nonumber
\end{align}
\end{widetext}
where all operators are in the Heisenberg representation. These expressions will be employed as an independent check on the correctness of the single-particle Green's function that we derive next.

\section{Time Dependence}
\label{sec:timedep}
\subsection{Time Evolution Operator}
To obtain the Green's functions of Eqs.~(\ref{gfdef}) and  (\ref{2partgf}), we must first determine the time evolution operator. To do this, we start by switching to the interaction representation with respect to the phonon part of the Hamiltonian and following the technique of Ref.~\onlinecite{jimHolstein}. In the interaction representation, we factorize the time evolution operator into two pieces via
\begin{equation}
    \hat{U}(t,t_{min})=e^{-i\hat{a}^\dagger\hat{a}\omega (t-t_{min})}\mathcal{T}_t\exp\left [-i\int_{t_{min}}^t\!\!\!dt^\prime\hat{\mathcal H}_I(t^\prime)\right ]
\end{equation}
where the interaction Hamiltonian has two pieces $\hat{\mathcal H}_I=\hat{\mathcal H}_I^{el-ph}+\hat{\mathcal H}_I^{coul}$. The electron-phonon coupled piece acquires an additional time dependence through the time dependence of the phonon raising and lowering operators
\begin{equation}
    \hat{\mathcal H}^{el-ph}_I(t)=\frac{g(t)}{\sqrt{2m\omega}}\big(e^{i\omega t} \hat{a}^{\dagger}+e^{-i\omega t}
    \hat{a}\big)\big(\hat{n}_{\sigma}+\hat{n}_{\bar{\sigma}}\big),
\end{equation}
while the Coulomb piece is unaltered ($\hat{\mathcal H}_I^{coul}(t)=U(t)\hat{n}_\uparrow\hat{n}_\downarrow$). These two pieces commute with each other, so the time evolution can be further factorized into a piece from the Coulomb part of the Hamiltonian and a piece from the electron-phonon part. The Coulomb piece is also simple, because it requires no time ordering, and is given by $\exp\left (-i\hat{n}_\uparrow\hat{n}_\downarrow\int_{t_{min}}^t\!\!\!dt^\prime U(t^\prime)\right  )$.
The electron-phonon piece is more complicated because the electron-phonon interaction Hamiltonian at two different times does not commute with itself:
\begin{equation}
    \left [\hat{\mathcal H}^{el-ph}_I(t),\hat{\mathcal H}^{el-ph}_I(t')\right ]=\frac{g(t)g(t^\prime)}{m\omega}i\sin\omega(t-t^\prime) (\hat{n}_{\uparrow}+\hat{n}_{\downarrow})^2,
\end{equation}
but the commutator does commute with $\hat{\mathcal H}^{el-ph}_I(t^{\prime\prime})$. This allows us to solve for the time-evolution operator following the strategy used by Landau and Lifshitz\cite{landau_lifshitz} and by Gottfried\cite{gottfried} in solving the driven simple harmonic oscillator.  We first form the function
\begin{equation}
    \hat{w}_I(t,t_{min})=\int_{t_{min}}^t\,dt'\ \hat{\mathcal H}^{el-ph}_I(t'),
\end{equation}
which satisfies the following equal-time commutator:
\begin{eqnarray}
    \left [\hat{\mathcal H}^{el-ph}_I(t),\hat{w}_I(t,t_{min})\right ]&=&i\frac{g(t)}{m\omega}\int_{t_{min}}^t\!\!\!dt^\prime g(t^\prime)\sin\omega(t-t^\prime)\nonumber\\
    &\times&(\hat{n}_{\uparrow}+\hat{n}_{\downarrow})^2.
\end{eqnarray}
This commutator commutes with all other operators in the different terms in the Hamiltonian.
When these two conditions hold, we have the following operator identity:
\begin{align}
          &e^{i\hat{w}_I(t,t_{min})}\left(i\frac{\partial}{\partial t}-\hat{\mathcal H}^{el-ph}_{I}(t)\right )e^{-i\hat{w}_I(t,t_{min})} \nonumber\\
    &~~~~~=i\frac{\partial}{\partial t}-\frac{i}{2}\left [\hat{w}_I(t,t_{min}),\hat{\mathcal H}^{el-ph}_{I}(t)\right ].
\end{align}
This allows us to evaluate the following time-ordered product directly as
\begin{align}
    &\mathcal{T}_t\exp\left[-i\int_{t_{min}}^t\!\!\!dt'\ \hat{\mathcal H}^{el-ph}_{I}(t')\right ]=\exp\left [-i\hat{w}_I(t,t_{min})\right ]\nonumber\\
&~~~~~\times\exp\left(\frac{1}{2}\int_{t_{min}}^t\!\!\!dt'\ \left [\hat{w}_I(t',t_{min}),\hat{\mathcal H}^{el-ph}_{I}(t')\right ]\right).
\end{align}
Applying this result to our specific Hamiltonian, yields the full time-evolution operator
\begin{align}
     &\hat{U}(t,t_{min})=\exp\left[-i\int_{t_{min}}^t\!\!\!dt'\left (-\mu \hat{n}_f+\omega \hat{a}^{\dagger}\hat{a}+\hat{n}_{\sig}\hat{n}_{\bar{\sigma}}U(t')\right )\right ]\nonumber\\ &\times\exp\left[-i\int_{t_{min}}^t\!\!\!dt' \frac{g(t')}{\sqrt{2m\omega}}(e^{i\omega t'}\hat{a}^{\dagger}+e^{-i\omega t'}\hat{a})\hat{n}_f\right]\\ &\times\exp\left[\frac{i}{2m\omega}\int_{t_{min}}^t\!\!\!dt' \int_{t_{min}}^{t'}\!\!\!dt''\ g(t')g(t'')\sin\omega(t'-t'')\hat{n}_f^2\right]
     \nonumber
\end{align}
where $\hat{n}_f=\hat{n}_{\uparrow}+\hat{n}_{\downarrow}$ is the total electron number operator. Note that because $\hat{n}_{\sigma}^2=\hat{n}_{\sigma}$, as required by the Pauli exclusion principle, we also have that $\hat{n}_f^2=(\hat{n}_{\uparrow}+\hat{n}_{\downarrow})^2=\hat{n}_{\uparrow}+\hat{n}_{\downarrow}+2\hat{n}_\uparrow \hat{n}_\downarrow$.
\subsection{Partition Function}
Calculating thermal averages, requires us to also determine the atomic partition function, which is defined to be
\begin{align}
    \mathcal{Z}&=\text{Tr}_{b,f}\big\{\hat{U}(t=-i\beta+t_{min},t_{min})\big\}\nonumber\\&=\text{Tr}_{b,f}\big\{e^{-\beta\hat{\mathcal H}(t_{min})}\big\},
\end{align}
where the trace is over product states composed of the tensor product of the Fermionic states $\ket{0},\ket{\uparrow},\ket{\downarrow},\ket{\uparrow\downarrow}$ with the Bosonic simple harmonic oscillator number states $|n\rangle=\frac{1}{\sqrt{n!}}\left (\hat{a}^\dagger\right )^n|0\rangle$. Evaluating the trace by using standard methods, we find
\begin{align}
\label{partition}
    \mathcal{Z}&=\Big(n(\omega)+1\Big)\left\{1+2e^{\beta\mu}\exp\left[\frac{\beta g^2(t_{min})}{2m\omega^2}\right]\right.\nonumber\\
    &+\left .\exp\left[\beta(2\mu-U(t_{min}))\right]\exp\left[\frac{2\beta g^2(t_{min})}{m\omega^2}\right]\right\}
\end{align}
with the Bose-Einstein distribution function given by
\begin{align}
\label{bose}
    n(\omega)=\frac{1}{e^{\beta\omega}-1}.
\end{align}
\subsection{Time Dependence of Field Operators}
Next we need the time dependence of the Fermionic field operators in the Heisenberg representation,
\begin{equation}
    \hat{c}_{\sigma}(t)=\hat{U}^{\dagger}(t,t_{min})\hat{c}_{\sigma}\hat{U}(t,t_{min})
\end{equation}
and similarly for the conjugate \begin{equation}
    \hat{c}_{\sigma}^{\dagger}(t)=\hat{U}^{\dagger}(t,t_{min})\hat{c}_{\sigma}^{\dagger}(t)\hat{U}(t,t_{min}).
\end{equation}
To evaluate these operator expressions, we note that the only operator in $\hat{U}(t,t_{min})$ that does not commute with $\hat{c}_{\sigma}$ is $\hat{n}_{\sigma}$. Hence, we need to compute
\begin{equation}
    f(\alpha)=e^{i\alpha \hat{n}_{\sigma}}\hat{c}_{\sigma}e^{-i\alpha \hat{n}_{\sigma}}.
\end{equation}
Differentiating gives us
\begin{equation}
    \frac{df(\alpha)}{d\alpha}=e^{i\alpha \hat{n}_{\sigma}}[\hat{n}_{\sigma},\hat{c}_{\sigma}]e^{-i\alpha \hat{n}_{\sigma}} =-if(\alpha),
\end{equation}
which then can be integrated to yield
\begin{equation}
    f(\alpha)=e^{-i\alpha}\hat{c}_{\sigma}.
\end{equation}\\
Hence, we find that\\
\begin{widetext}
\begin{align}
\hat{c}_{\sigma}(t)&=\exp\left[-i\int_{t_{min}}^t\!\!\!dt'\left (-\mu +\hat{n}_{\bar{\sigma}} U(t')\right )-\frac{i}{\sqrt{2m\omega}}\int_{t_{min}}^t\!\!\!dt'\ g(t')\left (e^{i\omega t'}\hat{a}^{\dagger}+e^{-i\omega t'}\hat{a}\right )\right]\nonumber\\&\times\exp\left[\frac{i}{2m\omega}\int_{t_{min}}^t\!\!\!dt'\ \int_{t_{min}}^{t'}\!\!\!dt''\ g(t')g(t'')\sin\omega(t'-t'')(1+2\hat{n}_{\bar{\sigma}})\right]\hat{c}_{\sigma}
\end{align}
and for the conjugate
\begin{align}
\hat{c}_{\sigma}^{\dagger}(t)&=\exp\left[i\int_{t_{min}}^t\!\!\!dt'\ \left (-\mu+\hat{n}_{\bar{\sigma}}U(t')\right ) +\frac{i}{\sqrt{2m\omega}}\int_{t_{min}}^t\!\!\!dt'\ g(t')\left (e^{i\omega t'}\hat{a}^{\dagger}+e^{-i\omega t'}\hat{a}\right )\right]\nonumber\\&\times\exp\left[-\frac{i}{2m\omega}\int_{t_{min}}^t\!\!\!dt'\ \int_{t_{min}}^{t'}\!\!\!dt''\ g(t')g(t'')\sin\omega(t'-t'')(1+2\hat{n}_{\bar{\sigma}})\right]\hat{c}_{\sigma}^{\dagger}.
\end{align}
\end{widetext}
With these time-dependent operators we are able to exactly evaluate the single- and two-particle Green's functions.
\section{Single-Particle Green's Function}
\label{sec:onepart}
To evaluate the single-particle Green's function in Eq.~(\ref{gfdef}), we must take a trace over the Hilbert space composed of a direct product of the Fermionic states with the simple harmonic oscillator states.
Newton's generalized binomial theorem enables us to sum the infinite degrees of freedom of the harmonic oscillator to obtain a closed form expression for the Bosonic trace. This again follows the technique used in Ref.~\onlinecite{jimHolstein} to evaluate the partition function. We outline this procedure briefly here, as it is used repeatedly in our derivations.  We have
\begin{align}
    g^R_{\sigma}(t_1,t_2)&=-\frac{i\theta(t_1-t_2)}{\mathcal Z}\nonumber\\
    &\times\text{Tr}_{b,f}\bigg\{e^{-\beta\hat{\mathcal H}(t_{min})}\{\hat{c}_{\sigma}(t_1),\hat{c}_{\sigma}^{\dagger}(t_2)\}_+\bigg\}
\end{align}
The Bosonic trace is calculated in the following way. First, we note that the single Bosonic operators always appear as summed pairs, which can be re-expressed by employing the Weyl form of the Baker-Campbell-Hausdorff identity as
\begin{equation}
e^{\alpha \hat{a}+\beta \hat{a}^{\dagger}}=e^{\alpha \hat{a}}e^{\beta \hat{a}^{\dagger}}e^{-\frac{1}{2}\alpha\beta[\hat{a},\hat{a}^{\dagger}]}.
\end{equation}
Using this expression, we can reduce the Bosonic trace to  constants (with respect to the Bosonic operators) multiplied by the form
\begin{equation}
    \text{Tr}_b\left\{e^{-\beta\omega \hat{a}^{\dagger}\hat{a}}e^{A\hat{a}}e^{B\hat{a}^{\dagger}}\right\},
\end{equation}
where $A$ and $B$ are constants \emph{with respect to the Bosonic operators}. Note that the simple harmonic oscillator states $\ket{n}$ are eigenvectors of the Bosonic number operator $\hat{n}=\hat{a}^{\dagger}\hat{a}$ with eigenvalue $n$, so we have
\begin{align}
     &\text{Tr}_b\left\{e^{-\beta\omega \hat{a}^{\dagger}\hat{a}}e^{A\hat{a}}e^{B\hat{a}^{\dagger}}\right\}=\sum_{n=0}^{\infty}e^{-\beta\omega  n}\langle n|e^{A\hat{a}}e^{B\hat{a}^{\dagger}}|n\rangle.
\end{align}
By expanding each exponential inside the expectation value in a power series, using the properties of the simple harmonic oscillator states and summing the resulting series with Newton's generalized binomial theorem, we evaluate the trace to be
\begin{align}
     &\text{Tr}_b\big\{e^{-\beta\omega \hat{a}^{\dagger}\hat{a}}e^{A\hat{a}}e^{B\hat{a}^{\dagger}}\big\}\nonumber\\
     &~~~~~=\frac{1}{1-e^{-\beta\omega}}\exp\left[\frac{1}{1-e^{-\beta\omega}}AB\right].
\end{align}
Using this result allows us to evaluate the Bosonic trace in the definition of the Green's function, which becomes
\begin{align}
\label{bosetrace}
    &g_{\sigma}^R(t_1,t_2)=\frac{-i\theta(t_1-t_2)\left(n(\omega)+1\right)}{\mathcal Z}e^{i\mu(t_1-t_2)}\nonumber\\ &~~~\times\exp\left[\left(-\frac{1}{2}-n(\omega)\right)\big|C(t_1)-C(t_2)\big|^2\right]\nonumber\\ &~~~\times
    \text{Tr}_f\bigg\{\exp\left[\frac{\beta g^2(t_{min})\hat{n}_f^2}{2m\omega^2}+\beta\mu \hat{n}_f-\beta U(t_{min})\hat{n}_{\uparrow}\hat{n}_{\downarrow}\right]\nonumber\\
    &~~~\times\exp\left[-i\int_{t_2}^{t_1}\,dt\ U(t)\hat{n}_{\bar{\sigma}}\right]\nonumber\\
    &~~~\times\exp\left[i(1+2\hat{n}_{\bar{\sigma}})\big(I(t_1)-I(t_2)\big)]\right]\nonumber\\
   &~~~\times \exp\left[\frac{2ig(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{C(t_1)-C(t_2)\big\}\hat{n}_f\right]\nonumber\\
    &~~~\times\left\{e^{i\text{Im}\{C^{*}(t_1)C(t_2)\}}\hat{c}_{\sigma}\hat{c}_{\sigma}^{\dagger}+ e^{-i\text{Im}\{C^{*}(t_1)C(t_2)\}}\hat{c}_{\sigma}^{\dagger}\hat{c}_{\sigma}\right\}\bigg\}
\end{align}
where
\begin{equation}
\label{Coft}
    C(t)=\frac{1}{\sqrt{2m\omega}}\int_{t_{min}}^t\!\!\!dt'\ g(t') e^{i\omega t'},
\end{equation}
and \begin{align}
\label{Ioft}
    I(t)=\frac{1}{2m\omega}\int_{t_{min}}^t\!\!\!dt' \int_{t_{min}}^{t'}\!\!\!dt''\ g(t')g(t'')\sin\omega(t'-t'')
\end{align}
are functions we define to make this result more readable.
We can now evaluate the Fermionic trace over the four possible states, with the result
\begin{widetext}
\begin{align}
\label{final1}
g_{\sigma}^R(t_1,t_2)&=-\frac{i\theta(t_1-t_2)\big (n(\omega)+1\big)}{\mathcal Z}e^{i\mu(t_1-t_2)}
\exp\bigg[\left(-\frac{1}{2}-n(\omega)\right)\big|C(t_1)-C(t_2)\big|^2\bigg]\nonumber\\&\times
\bigg\{e^{i\text{Im}\{C^{*}(t_1)C(t_2)\}}\bigg(\exp\big[i\big(I(t_1)-I(t_2)\big)\big]\nonumber\\
&+\exp\left[\beta(\mu+\frac{g^2(t_{min})}{2m\omega^2})\right] \exp\left[-i\int_{t_2}^{t_1}\,dt'\ U(t')\right]\exp\big[3i(I(t_1)-I(t_2))\big]
\exp\left[\frac{2ig(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{C(t_1)-C(t_2)\}\right ]\bigg)\nonumber\\
&+e^{-i\text{Im}\{C^{*}(t_1)C(t_2)\}}\bigg( \exp\left[\beta(\frac{\mu+g^2(t_{min})}{2m\omega^2})\right] \exp\big[i\big(I(t_1)-I(t_2)\big)\big] \exp\left[\frac{2ig(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{C(t_1)-C(t_2)\}\right]\nonumber\\
&+\exp\left[\frac{2\beta g^2(t_{min})}{m\omega^2}+\beta(2\mu-U(t_{min}))\right]\exp\left[-i\int_{t_2}^{t_1}\,dt\ U(t)\right]\exp\big[3i(I(t_1)-I(t_2))\big]\nonumber\\
&\times\exp\left[\frac{4ig(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{C(t_1)-C(t_2)\}\right]\bigg)\bigg\}
\end{align}
\end{widetext}

This result is quite complicated. But note that it is an exact expression, so it can be employed to describe the response of Holstein-Hubbard systems with arbitrary time-dependent interactions. It is remarkable that such a complex system has such a compact exact expression for the nonequilibrium Green's functions. We reiterate again that there are two keys to allowing this to occur: first the number operators for each spin are conserved and do not depend on time and second, the commutator of the electron-phonon interaction with itself at two different times is an operator that commutes with all remaining operators in the system. Nevertheless, the final result is not simple. It represents the complex dynamics of electron-phonon coupled systems in the atomic limit. In contrast, the atomic limit Green's functions for the Hubbard model are quite simple.

While it would be nice to have an independent check to verify that these results are correct, we are not aware of any way to do so for the full functions. But we do demonstrate below that this Green's function does satisfy the first four spectral moment sum rules. The first three are evaluated analytically by taking derivatives of Eq.~(\ref{bosetrace}), before we take the Fermionic trace, as in Eq.~(\ref{derivs}). The third moment is verified by taking numerical derivatives of the equation Eq.~(\ref{final1}). We also numerically check the first three moments to ensure we have taken the trace correctly in going from Eq.~(\ref{bosetrace}) to Eq.~(\ref{final1}). While this does not test everything, it is a quite stringent test of accuracy of the final results.
\section{Verification of Moments}
\label{sec:moms2}
\begin{figure}[b]
    \centering
    \includegraphics[width=0.4\textwidth,angle=270]{thirdmomHH.jpg}
    \caption{Comparison of the exact third moment  [given by Eq.~(\ref{mom3})] and the numerical third derivative of Eq.~(\ref{final1}) for non-constant electron-phonon and electron-electron coupling parameters {$g(t)=1-0.05e^{-0.01t}$} and {$U(t)=1-0.01t^2$}.}
    \label{fig:Third}
\end{figure}
Using the definition of the Green's function and Eq.~(\ref{derivs}), allows us to write the moments as derivatives of the causal Green's function with respect to relative time  \begin{equation}
    \mu_{\sigma}^{Rn}(t_{ave})=-\text{Im}\big\{i^n\frac{\partial^n}{\partial t_{rel}^n}g^R_{\sigma}(t_{ave}+\frac{t_{rel}}{2},t_{ave}-\frac{t_{rel}}{2})\big\}\big|_{t_{rel\to 0}}.
\end{equation}
Using this form, we easily find from Eq.~\ref{final1} that $\mu_{\sigma}^{R0}(t_{ave})=1$. By evaluating traces similar to those we evaluated to determine the partition function and the Green's function, we also find exact expressions for the time-dependent expectation values, which appear in the moments of Eqs.~(\ref{mom1}-\ref{mom3}). Following Ref.~\onlinecite{jimHolstein}, we have following two expectation values:
\begin{widetext}
    \begin{equation}
    \langle \hat{x}(t)\rangle =-\langle \hat{n}_f\rangle \left (\frac{g(t_{min})}{m\omega^2}\cos\big(\omega (t-t_{min})\big)+\text{Re}\left\{ie^{-i\omega t}\int_{t_{min}}^t\!\!\!dt'\ \frac{g(t')}{m\omega}e^{i\omega t'}\right\}\right)
\end{equation}
and \begin{equation}
\label{nf}
    \langle \hat{n}_f\rangle =\frac{2e^{\beta\mu}\exp\left[\frac{\beta g^2(t_{min})}{2m\omega^2}\right]+2e^{\beta(2\mu-U(t_{min}))}\exp\left[\frac{2\beta g^2(t_{min})}{m\omega^2}\right]}{1+2e^{\beta\mu}\exp\left[\frac{\beta g^2(t_{min})}{2m\omega^2}\right]
    +\exp\left[\beta(2\mu-U(t_{min}))\right]\exp\left[\frac{2\beta g^2(t_{min})}{m\omega^2}\right]}
\end{equation}
\end{widetext}
The first is time dependent, while the second is time independent (because electron number is conserved in the model).
The remaining expectation values needed to evaluate the analytic expressions for the moments are given in Appendix~\ref{sec:a}. By taking derivatives of Eq.~(\ref{bosetrace}) and simplifying we find \begin{equation}
    \mu_{\sigma}^{R1}=-\mu+U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle +g(t_{ave})\langle \hat{x}(t_{ave})\rangle,
\end{equation}
in agreement with Eq.~(\ref{mom1})) and
\begin{widetext}
\begin{align}
    \mu_{\sigma}^{R2}(t_{ave})&=\big(\mu-U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle -g(t_{ave})\langle \hat{x}(t_{ave})\rangle\big)^2+U(t_{ave})^2
    \big[\langle \hat{n}_{\bar{\sigma}}\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^2\big]\nonumber\\
    &+g^2(t_{ave})\big[\langle \hat{x}^2(t_{ave})\rangle-\langle \hat{x}(t_{ave})\rangle^2\big]+2U(t_{ave})g(t_{ave})\big[\langle \hat{n}_{\bar{\sigma}}\hat{x}(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle \langle \hat{x}(t_{ave})\rangle\big],
\end{align}
\end{widetext}
in agreement with Eq.~(\ref{mom2}), proving that the first three moments of the single-particle Green's function hold exactly. For the third moment, taking derivatives by hand becomes too cumbersome, so we resort to taking numerical derivatives to verify that the third moment holds. In Fig.~\ref{fig:Third}, we plot this comparison in the case where $t_{min}=0$, $U(t)=1-0.01t^2$ and $g(t)=1-0.05\exp[-0.1t]$. One can clearly see the third moment is verified numerically. While we plot just one example, this agreement held for all combinations of parameters that we tested, providing strong evidence that the third-moment sum rule is also satisfied by the single-particle Green's function in Eq.~(\ref{final1}).
While this exercise is not definitive, it is hard to imagine how one can satisfy all of these moment sum rules and still have an error (of course, it could still occur with the higher-order terms, but trying to evaluate higher-order moments is also too cumbersome).

\section{Two Particle Green's Function}
\label{sec:twopart}
For the two particle Green's function, there are no known sum rules to compare our results against, so we no longer have an independent check. The definition of the two-particle Green's function is given in Eq.~(\ref{2partgf}). There are two possibilities for the spin $\sigma'$, $\sigma'=\sigma$ or $\sigma'=\bar{\sigma}$.
\subsection{Same Spin Case}
For the case where the spins are the same, we find a number of terms in the time ordering vanish, because $\hat{c}^2_{\sigma}=\hat{c}_{\sigma}^{\dagger 2}=0$ due to the Pauli exclusion principle. In this case, only 8 of the 24 possible time-orderings are non-vanishing. When we expand the time-ordering operator for the same spin Green's function, we obtain
\begin{widetext}
\begin{align}
\label{samespin}
    \mathcal{G}_{\sigma\sigma}(t_0,t_1,t_2,t_3)&=-\langle {\mathcal T}_t\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}(t_1)\hat{c}^{\dagger}_{\sigma}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle
    \\&=\theta(t_0,t_2,t_1,t_3)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle-\theta(t_0,t_3,t_1,t_2)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\rangle \nonumber\\
    &+\theta(t_1,t_3,t_0,t_2)\langle \hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\rangle -\theta(t_1,t_2,t_0,t_3)\langle \hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle \nonumber\\
    &+\theta(t_2,t_0,t_3,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\rangle -\theta(t_2,t_1,t_3,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\rangle \nonumber\\
    &+\theta(t_3,t_1,t_2,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\rangle -\theta(t_3,t_0,t_2,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\rangle. \nonumber
\end{align}

The expectation values are evaluated in the same fashion as we did for the single-particle case. For example, the first term, which is nonzero when $t_0>t_2>t_1>t_3$, is given by
\begin{align}
    &\langle \hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\sigma}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle
    =
    \frac{\big(n(\omega)+1\big)}{\mathcal Z}e^{-i\mu(-t_0+t_2-t_1+t_3)}
    \exp\left[\left(-\frac{1}{2}-n(\omega)\right)\big|C(t_2)-C(t_0)+C(t_1)-C(t_3)\big|^2\right]\\ &~~~~~~\times\exp\big[i\big(I(t_0)-I(t_2)+I(t_1)-I(t_3)\big)\big]
    \exp[P(-t_0,t_2,-t_1,t_3)]\nonumber\\
    &~~~~~~\times\bigg\{1+\exp\left[\beta\left(\mu+\frac{g^2(t_{min})}{2m\omega^2}\right)\right] \exp\left[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\}\right]\nonumber\\ &~~~~~~\times\exp\left[-i\left(-2I(t_0)+\int_{t_{min}}^{t_0}\!\!\!dt' U(t')+2I(t_2)-\int_{t_{min}}^{t_2}\!\!\!dt'U(t')
    -2I(t_1)+\int_{t_{min}}^{t_1}\!\!\!dt' U(t')+2I(t_3)-\int_{t_{min}}^{t_3}\!\!\!dt' U(t')\right)\right]\bigg\}\nonumber
\end{align}
where $C(t)$ and $I(t)$ are defined in Eqs.~(\ref{Coft}) and (\ref{Ioft}) and $P(t,t', \bar{t},\bar{t}')$ is a prefactor, which results from the repeated application of the Baker-Campbell-Hausdorff theorem, and is defined by the following expression
\begin{align}
        &P((\pm)_1t,(\pm)_2t',(\pm)_3\bar{t},(\pm)_4\bar{t}')=
        \exp\bigg[-i\text{Im}\big\{(\pm)_1(\pm)_2C^{*}(t)C(t')
        +(\pm)_3(\pm)_4C^{*}(\bar{t})C(\bar{t}')\\&~~~~~~+(\pm)_1(\pm)_3C^{*}(t)C(\bar{t})
        +(\pm)_1(\pm)_4C^{*}(t)C(\bar{t}')+(\pm)_2(\pm)_3C^{*}(t')C(\bar{t})
        +(\pm)_2(\pm)_4C^{*}(t')C(\bar{t}')\big\}\bigg]\nonumber
\end{align}
\end{widetext}
where the $(\pm)_i$ is a $+$ if the time $i$ is associated with a $c^{\dagger}$ operator, $-$ for a $c$ operator.
The expressions for other time orderings are given in appendix~\ref{sec:b}.
\subsection{Different Spin Case}
In the case where $\sigma'=\bar{\sigma}$, none of the terms in the time-ordering vanish. Then the Green's function for the mixed spin case is given by
\begin{widetext}
\begin{align}
\label{differentspin}
\mathcal{G}_{\sigma\bar{\sigma}}(t_0,t_1,t_2,t_3)&=
-\theta(t_0,t_1,t_2,t_3)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle+\theta(t_1,t_0,t_2,t_3)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle \nonumber\\
&-\theta(t_1,t_0,t_3,t_2)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle +\theta(t_0,t_1,t_3,t_2)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\rangle \nonumber\\
&-\theta(t_0,t_2,t_3,t_1)\langle \hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle +\theta(t_0,t_2,t_1,t_3)\langle \hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\rangle \nonumber\\
&-\theta(t_0,t_3,t_1,t_2)\langle \hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle +\theta(t_0,t_3,t_2,t_1)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle \nonumber\\
&-\theta(t_1,t_3,t_2,t_0)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}(t_0)\rangle +\theta(t_1,t_2,t_3,t_0)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\rangle \nonumber\\
&-\theta(t_1,t_2,t_0,t_3)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle +\theta(t_1,t_3,t_0,t_2)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\rangle \nonumber\\
&-\theta(t_2,t_0,t_1,t_3)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle +\theta(t_2,t_1,t_0,t_3)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle \nonumber\\
&-\theta(t_2,t_1,t_3,t_0)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\rangle +\theta(t_2,t_3,t_1,t_0)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\rangle \nonumber\\
&-\theta(t_2,t_3,t_0,t_1)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\rangle +\theta(t_2,t_0,t_3,t_1)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle \nonumber\\
&-\theta(t_3,t_1,t_0,t_2)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle +\theta(t_3,t_0,t_1,t_2)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle \nonumber\\
&-\theta(t_3,t_2,t_1,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\rangle +\theta(t_3,t_2,t_0,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\rangle \nonumber\\
&-\theta(t_3,t_0,t_2,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\rangle +\theta(t_3,t_1,t_2,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\rangle.
\end{align}
For example, the term with $t_0>t_2>t_1>t_3$ is given by
\begin{align}
           &\langle c_{\sigma}(t_0)c_{\bar{\sigma}}^{\dagger}(t_2)c_{\bar{\sigma}}(t_1)c_{\sigma}^{\dagger}(t_3)\rangle = \frac{\big(n(\omega)+1\big)}{\mathcal Z}e^{-i\mu(-t_0+t_2-t_1+t_3)}
    \exp\left[\left(-\frac{1}{2}-n(\omega)\right)\big|C(t_2)-C(t_0)+C(t_1)-C(t_3)\big|^2\right]\nonumber\\
    &~~~~~\times\exp\big[i\big(I(t_0)-I(t_2)+I(t_1)-I(t_3)\big)\big]
           e^{P(-t_0,t_2,-t_1,t_3)}
           \exp\left[\frac{\beta g^2(t_{min})}{2m\omega^2}\right]\exp\big[\beta\mu\big]
           \nonumber\\ &~~~~~\times\exp\left[-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(-t_0)+C(t_2)-C(t_1)+C(t_3)\}\right]\nonumber\\
      &~~~~~\times\exp\bigg[-i\big(-2I(t_0)+\int_{t_{min}}^{t_0}\!\!\!dt' U(t')+2I(t_2)-\int_{t_{min}}^{t_2}\!\!\!dt' U(t')
    -2I(t_1)+\int_{t_{min}}^{t_1}\!\!\!dt' U(t')+2I(t_3)-\int_{t_{min}}^{t_3}\!\!\!dt' U(t')\big)\bigg]\bigg\}
       \end{align}
      \end{widetext}
While there are many terms, the evaluation is again straightforward, but lengthy. Results for the other time orderings in the mixed spin case are also presented in Appendix~\ref{sec:b}.
\section{Conclusion}
\label{sec:conclusion}
In this work, we have found an exact result for the single-particle and two-particle correlation functions of the nonequilibrium Holstein-Hubbard model (with time dependent interactions) in the atomic limit. This includes all of the complex dynamical effects associated with the coupling of electrons to quantum phonons and the possible change in time of the electron-electron and electron-phonon interactions. Such dynamics are, of course, quite complex, so the results are cumbersome. While there is no independent check on their correctness, we did partially verify them by comparing them to exact results of the first four nonequilibrium spectral moments of the single-particle retarded Green's function, finding exact agreement.  We also calculated the two-particle Green's functions for this model. Both of these Green's functions are needed to calculate the lattice Green's functions within a strong-coupling-expansion framework. We believe such future work will eventually  provide insight into the behavior of complex materials in pump-probe spectroscopy experiments, especially those that have bipolarons. In the future, it may be interesting to also expand these results to the case with a time-dependent  phonon frequency and to consider non-linear couplings of the electrons to the phonons. This would be even more impactful if one can properly determine how the change in electron occupancy during a pump modifies the interaction parameters self-consistently. Perhaps these exact results could be employed within a more complicated formalism to do just that for systems that are near the atomic limit. But, even as standalone results, these expressions provide potential benchmarks to nonequilibrium many-body calculations that examine the behavior of strongly coupled electron-phonon systems near the atomic limit.

\section*{Acknowledgements} \label{sec:acknowledgements}
This work was supported by
the Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Contract No.~DE-SC0019126.
J.~K.~F. was also supported by the McDevitt bequest at Georgetown. We acknowledge useful discussions with K.~Najafi.

\section*{Acknowledgements} \label{sec:acknowledgements}
This work was supported by
the Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Contract No.~DE-SC0019126.
J.~K.~F. was also supported by the McDevitt bequest at Georgetown. We acknowledge useful discussions with K.~Najafi.
\section{Expectation Values}
    \label{sec:a}
    In this Appendix, we summarize the calculation of the trace needed for different time-dependent expectation values.
This includes the terms $\langle \hat{n}_{\bar{\sigma}}\rangle, \langle \hat{n}_{\sigma}\hat{n}_{\bar{\sigma}}\rangle,$ $(\langle \hat{x}^2(t)\rangle-\langle \hat{x}(t)\rangle^2), (\langle \hat{x}^3(t)\rangle-\langle \hat{x}(t)\rangle ^3$, $\langle \hat{n}_{\bar{\sigma}}\hat{x}(t)\rangle$ and $\langle \hat{n}_{\bar{\sigma}}\hat{x}^2(t)\rangle$. They are evaluated by using Newton's generalized binomial theorem and properties of the simple-harmonic-oscillator states, as outlined in the main paper. To begin,
we have \begin{equation}
    \langle \hat{n}_{\bar{\sigma}}\rangle =\frac{e^{\beta\mu}\exp\left[\frac{g^2(t_{min})\beta}{2m\omega^2}\right]+e^{\beta\left(2\mu-U\left(t_{min}\right)\right)}\exp\left[\frac{2g^2(t_{min})\beta}{m\omega^2}\right]}{1+2e^{\beta\mu}\exp\left[\frac{g^2(t_{min})\beta}{2m\omega^2}\right]
    +e^{\beta\left(2\mu-U\left(t_{min}\right)\right)}\exp \left[\frac{2g^2(t_{min})\beta}{m\omega^2}\right]}
\end{equation}
so we see, by comparison with Eq.~(\ref{nf}), that we have $\langle \hat{n}_{\bar{\sigma}}\rangle=\frac{1}{2}\langle \hat{n}_{\sigma}+\hat{n}_{\bar{\sigma}}\rangle=\frac{1}{2}\langle\hat{n}_f\rangle$. In addition,
\begin{widetext}\begin{align}
    \langle \hat{n}_{\sigma}\hat{n}_{\bar{\sigma}}\rangle =\frac{e^{\beta\left(2\mu-U\left(t_{min}\right)\right)}\exp\left[\frac{2g^2(t_{min})\beta}{m\omega^2}\right]}{1+2e^{\beta\mu}\exp\left[\frac{g^2(t_{min})\beta}{2m\omega^2}\right]
    +\exp\left[\beta\left(2\mu-U(t_{min})\right)\right]\exp\left[\frac{2g^2(t_{min})\beta}{m\omega^2}\right]},
\end{align}
so we have \begin{align}
    \begin{split}
    \langle \hat{n}_{\bar{\sigma}}\hat{x}(t)\rangle = -\left(\frac{g(t_{min})}{m\omega^2}\cos\left(\omega\left( t-t_{min}\right)\right)+\text{Re}\left\{ie^{-i\omega t}\int_{t_{min}}^t\!dt'\ \frac{g(t')}{m\omega}e^{i\omega t'}\right\}\right)\left(\langle \hat{n}_{\sigma}\hat{n}_{\bar{\sigma}}\rangle +\langle \hat{n}_{\bar{\sigma}}\rangle\right).
\end{split}
\end{align}
Next, \begin{equation}
    \langle \hat{n}_{\bar{\sigma}}\hat{x}^2(t)\rangle=\left(-\frac{g(t_{min})}{m\omega^2}\cos\left(\omega\left( t-t_{min}\right)\right)-\text{Re}\left\{ie^{-i\omega t}\int_{t_{min}}^t\!dt'\ \frac{g(t')}{m\omega}e^{i\omega t'}\right\}\right)^2\left(\langle \hat{n}_{\bar{\sigma}}\rangle+3\langle \hat{n}_{\sigma}\hat{n}_{\bar{\sigma}}\rangle\right)+\frac{1}{2m\omega}\coth(\frac{\beta\omega}{2})\langle \hat{n}_{\bar{\sigma}}\rangle
\end{equation}

In addition, \begin{align}
    \langle \hat{x}^2(t)\rangle -\langle \hat{x}(t)\rangle^2&=
\left(\frac{g(t_{min})}{m\omega^2}\cos\left(\omega\left( t-t_{min}\right)\right)+\text{Re}\left\{ie^{-i\omega t}\int_{t_{min}}^t\!dt'\ \frac{g(t')}{m\omega}e^{i\omega t'}\right\}\right)^2\\
&\times\left(\langle \hat{n}_{\uparrow}+\hat{n}_{\downarrow}\rangle -\langle \hat{n}_{\uparrow}+\hat{n}_{\downarrow}\rangle^2 +2\langle \hat{n}_{\uparrow}\hat{n}_{\downarrow}\rangle\right)+\frac{1}{2m\omega}\coth\left(\frac{\beta\omega}{2}\right)
\end{align}
And finally,  \begin{align} \begin{split}
\langle \hat{x}^3(t)\rangle -\langle \hat{x}(t)\rangle^3 &= -\left(\frac{g(t_{min})}{m\omega^2}\cos\left(\omega\left( t-t_{min}\right)\right)+\text{Re}\left\{ie^{-i\omega t}\int_{t_{min}}^t\!dt'\ \frac{g(t')}{m\omega}e^{i\omega t'}\right\}\right)^3\left(\langle \hat{n}_{\uparrow}+\hat{n}_{\downarrow}\rangle -\langle \hat{n}_{\uparrow}+\hat{n}_{\downarrow}\rangle^3+6\langle \hat{n}_{\uparrow}\hat{n}_{\downarrow}\rangle \right)\\&+\frac{3}{2m\omega}\coth\left(\frac{\beta\omega}{2}\right) \langle \hat{x}(t)\rangle
\end{split}
\end{align}
\end{widetext}
These are all the expectation values required to evaluate the spectral moments up to third order.

\section{Two-particle Green's function}
\label{sec:b}
Here we will enumerate all the terms which appear in the two-particle Green's function, defined by
\begin{equation}
    \mathcal{G}_{\sigma\sigma'}(t_0,t_1,t_2,t_3)=-\langle {\mathcal T}_t \hat{c}_{\sigma}(t_1)\hat{c}_{\sigma'}(t_1)\hat{c}_{\sigma'}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle
\end{equation}
To begin, there is a time-dependent factor which is the same for every term, so we'll define it here.
\begin{align}
   & A(t_0,t_1,t_2,t_3)\equiv A =\frac{\big (n(\omega)+1\big )}{Z}e^{i\mu(t_0+t_1-t_2-t_3)}\nonumber\\ &\times\exp\left[\left(-\frac{1}{2}-n(\omega)\right)\left|C(t_2)+C(t_3)-C(t_0)-C(t_1)\right|^2\right]\nonumber\\ &\times\exp\big[i\big(I(t_0)+I(t_1)-I(t_2)-I(t_3)\big)\big]
\end{align}
where again $n(\omega), C(t)$ and $I(t)$ are given by Eqs. ~(\ref{bose}), ~(\ref{Coft}), and ~(\ref{Ioft}) respectively.
Next define
\begin{align}
    F(t)\equiv2I(t)-\int_{t_{min}}^t\!dt'\ U(t')
\end{align}
and the prefactor, which depends on the four times (and whether each time is associated with a creation or an annihilation operator),
\begin{widetext}
\begin{align}
        P((\pm)_1t,(\pm)_2t',(\pm)_3\bar{t},(\pm)_4\bar{t}') & = \exp\big[-i\text{Im}\{(\pm)_1(\pm)_2C^{*}(t)C(t')
        +(\pm)_3(\pm)_4C^{*}(\bar{t})C(\bar{t}')+(\pm)_1(\pm)_3C^{*}(t)C(\bar{t})\nonumber
       \\& +(\pm)_1(\pm)_4C^{*}(t)C(\bar{t}')+(\pm)_2(\pm)_3C^{*}(t')C(\bar{t})
        +(\pm)_2(\pm)_4C^{*}(t')C(\bar{t}')\}\big],
\end{align}
\end{widetext}
where the $+$ sign is for a time associated with a creation operator, and the $-$ sign is for a time associated with an annihilation operator. Finally, we note that the sign of the term is negative (positive) if the parity of the permutation which brings the times back to the order $t_0>t_1>t_2>t_3$ is even (odd).
\subsection{Same Spin Case}
For the case where $\sigma'=\sigma$, only time orderings which do not place two creation or annihilation operators next to one another survive. These are the following eight time orderings:  \begin{enumerate}
    \item $t_0>t_2>t_1>t_3$
    \item $t_0>t_3>t_1>t_2$
    \item $t_1>t_2>t_0>t_3$
    \item $t_1>t_3>t_0>t_2$
    \item $ t_2>t_1>t_3>t_0$
    \item $t_2>t_0>t_3>t_1$
    \item $t_3>t_1>t_2>t_0$
    \item $t_3>t_0>t_2>t_1$
\end{enumerate}
The same-spin two-particle Green's function is given by Eq.~(\ref{samespin}):
\begin{widetext}
\begin{align}
\begin{split}
    \mathcal{G}_{\sigma\sigma}(t_0,t_1,t_2,t_3)&=-\langle {\mathcal T}_t\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}(t_1)\hat{c}^{\dagger}_{\sigma}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle
    \\&=\theta(t_0,t_2,t_1,t_3)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle-\theta(t_0,t_3,t_1,t_2)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\rangle \nonumber\\
    &+\theta(t_1,t_3,t_0,t_2)\langle \hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\rangle -\theta(t_1,t_2,t_0,t_3)\langle \hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle \nonumber\\
    &+\theta(t_2,t_0,t_3,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\rangle -\theta(t_2,t_1,t_3,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\rangle \nonumber\\
    &+\theta(t_3,t_1,t_2,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\rangle -\theta(t_3,t_0,t_2,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\rangle. \nonumber
    \end{split}
\end{align}
The expectation values are as follows:
    \begin{enumerate}
    \item $t_0>t_2>t_1>t_3$: \begin{align}
       & \langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle = A\exp\left[P(-t_0,t_2,-t_1,t_3)\right] \bigg\{1+\exp\left[\beta\left(\mu+\frac{g^2(t_{min})}{2m\omega^2}\right)\right] \nonumber \\
       & \times\exp\left[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\left\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\right\}\right]\exp\left[-i\left(-F(t_0)+F(t_2)-F(t_1)-F(t_3)\right)\right]\bigg\}\nonumber
    \end{align}
    \item $t_0>t_3>t_1>t_2$:\begin{align}\begin{split}
         &- \langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle = -A\exp\big[P(-t_0,t_3,-t_1,t_2)\big]\big\{1+\exp\big[\beta\big(\mu+\frac{g^2(t_{min})}{2m\omega^2}\big)\big]\\
        &\times\exp\big[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]\exp\big[-i\big(-F(t_0)+F(t_2)-F(t_1)-F(t_3)\big)\big]\big\}
    \end{split}
    \end{align}
    \item $t_1>t_2>t_0>t_3$:
    \begin{align}\begin{split}
        & - \langle \hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle = -A\exp\big[P(-t_1,t_2,-t_0,t_3)\big]\big\{1+\exp\big[\beta\big(\mu+\frac{g^2(t_{min})}{2m\omega^2}\big)\big]\\
        &\times\exp\big[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]\exp\big[-i\big(-F(t_0)+F(t_2)-F(t_1)-F(t_3)\big)\big]\big\}
    \end{split}
    \end{align}
    \item $t_1>t_3>t_0>t_2$: \begin{align}\begin{split}
          \langle &\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle = A\exp\big[P(-t_1,t_3,-t_0,t_2)\big]\big\{1+\exp\big[\beta\big(\mu+\frac{g^2(t_{min})}{2m\omega^2}\big)\big]\\
        &\times\exp\big[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]\exp\big[-i\big(-F(t_0)+F(t_2)-F(t_1)-F(t_3)\big)\big]\big\}
    \end{split}
    \end{align}
    Note that for these orderings all that changed is the ordering of the times in the prefactor $P$. We see something similar for the next four orderings.
    \item $t_2>t_1>t_3>t_0:$ \begin{align}
        \begin{split}
            &-\langle \hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\rangle = -A\exp\big[P(t_2,-t_1,t_3,-t_0)\big]\bigg\{\exp\big[\beta\big(\mu+\frac{g^2(t_{min})}{2m\omega^2}\big)\big]\\
           & \times\exp\big[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]+\exp\big[\beta\big(2\mu-U(t_{min})+\frac{2g^2(t_{min})}{m\omega^2}\big)\big]\\
            &\times\exp\big[-i\big(-F(t_0)+F(t_2)-F(t_1)+F(t_2)\big)\big]\exp\big[-i\frac{4g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]\bigg\}
        \end{split}
    \end{align}
    \item $t_2>t_0>t_3>t_1$:
    \begin{align}\begin{split}
            &\langle\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\rangle = A\exp\big[P(t_2,-t_0,t_3,-t_1)\big]\bigg\{\exp\big[\beta\big(\mu+\frac{g^2(t_{min})}{2m\omega^2}\big)\big]\\
            &\times\exp\big[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]+\exp\big[\beta\big(2\mu-U(t_{min})+\frac{2g^2(t_{min})}{m\omega^2}\big)\big]\\
            &\times\exp\big[-i\big(-F(t_0)+F(t_2)-F(t_1)+F(t_2)\big)\big]\exp\big[-i\frac{4g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]\bigg\}
    \end{split}\end{align}
    \item $t_3>t_1>t_2>t_0$:\begin{align}
        \begin{split}
            &\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\rangle = A\exp\big[P(t_3,-t_1,t_2,-t_0)\big]\bigg\{\exp\big[\beta\big(\mu+\frac{g^2(t_{min})}{2m\omega^2}\big)\big]\\
            &\times\exp\big[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]+\exp\big[\beta\big(2\mu-U(t_{min})+\frac{2g^2(t_{min})}{m\omega^2}\big)\big]\\
            &\times\exp\big[-i\big(-F(t_0)+F(t_2)-F(t_1)+F(t_2)\big)\big]\exp\big[-i\frac{4g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]\bigg\}
        \end{split}
    \end{align}
    \item $t_3>t_0>t_2>t_1$:\begin{align}
        \begin{split}
           &- \langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\rangle = -A\exp\big[P(t_3,-t_0,t_2,-t_1)\big]\bigg\{\exp\big[\beta\big(\mu+\frac{g^2(t_{min})}{2m\omega^2}\big)\big]\\
           & \times\exp\big[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]+\exp\big[\beta\big(2\mu-U(t_{min})+\frac{2g^2(t_{min})}{m\omega^2}\big)\big]\\
            &\times\exp\big[-i\big(-F(t_0)+F(t_2)-F(t_1)+F(t_2)\big)\big]\exp\big[-i\frac{4g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]\bigg\}
        \end{split}
    \end{align}
\end{enumerate}
So for this case, unlike the mixed spin case considered next, we can write a reasonably concise expression for the two particle Green's function:
\begin{align}
\begin{split}
    &\mathcal{G}_{\sigma\sigma}(t_0,t_1,t_2,t_3)=A\bigg(\bigg\{1+\exp\big[\beta\big(\mu+\frac{g^2(t_{min})}{2m\omega^2}\big)\big]\exp\big[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]\\&\times\exp\big[-i\big(-F(t_0)+F(t_2)-F(t_1)-F(t_3)\big)\big]\bigg\}\\
    &\times\bigg\{\theta(t_0,t_2,t_1,t_3)\exp\big[P(-t_0,t_2,-t_1,t_3)\big]-\theta(t_0,t_3,t_1,t_2)\exp\big[P(-t_0,t_3,-t_1,t_2)\big]\\&+\theta(t_1,t_3,t_0,t_2)\exp\big[P(-t_1,t_3,-t_0,t_2)\big]-\theta(t_1,t_2,t_0,t_3)\exp\big[P(-t_1,t_2,-t_0,t_3)\big]\bigg\}\\
    &+\bigg\{\exp\big[\beta\big(\mu+\frac{g^2(t_{min})}{2m\omega^2}\big)\big]\exp\big[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]+\exp\big[\beta\big(2\mu-U(t_{min})+\frac{2g^2(t_{min})}{m\omega^2}\big)\big]\\
            &\times\exp\big[-i\big(-F(t_0)+F(t_2)-F(t_1)+F(t_2)\big)\big]\exp\big[-i\frac{4g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\big\}\big]\bigg\}\\& \times\bigg\{\theta(t_2,t_0,t_3,t_1)\exp\big[P(t_2,-t_0,t_3,-t_1)\big]-\theta(t_2,t_1,t_3,t_0)\exp\big[P(t_2,-t_1,t_3,-t_0)\big]\\&+\theta(t_3,t_1,t_2,t_0)\exp\big[P(t_3,-t_1,t_2,-t_0)\big]-\theta(t_3,t_0,t_2,t_1)\exp\big[P(t_3,-t_0,t_2,-t_1)\big]\bigg\}\bigg)\end{split}
\end{align}
\end{widetext}
where all the functions which appear are defined above.
\subsection{Mixed Spin Case}
For the case when $\sigma'=\bar{\sigma}$, none of the orderings vanish and each of the 24 terms has a slightly different structure. We will enumerate each of the terms here. To begin the mixed spin Green's function is defined as\begin{widetext} \begin{align}
\begin{split}
\mathcal{G}_{\sigma\bar{\sigma}}(t_0,t_1,t_2,t_3)=
-\theta(t_0,t_1,t_2,t_3)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)c^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle+\theta(t_1,t_0,t_2,t_3)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle \\
-\theta(t_1,t_0,t_3,t_2)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle +\theta(t_0,t_1,t_3,t_2)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)c^{\dagger}_{\sigma}(t_3)c^{\dagger}_{\bar{\sigma}}(t_2)\rangle \\
-\theta(t_0,t_2,t_3,t_1)\langle \hat{c}_{\sigma}(t_0)c^{\dagger}_{\bar{\sigma}}(t_2)c^{\dagger}_{\sigma}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle +\theta(t_0,t_2,t_1,t_3)\langle \hat{c}_{\sigma}(t_0)c^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\bar{\sigma}}(t_1)c^{\dagger}_{\sigma}(t_3)\rangle \\
-\theta(t_0,t_3,t_1,t_2)\langle \hat{c}_{\sigma}(t_0)c^{\dagger}_{\sigma}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle +\theta(t_0,t_3,t_2,t_1)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle \\
-\theta(t_1,t_3,t_2,t_0)\langle \hat{c}_{\bar{\sigma}}(t_1)c^{\dagger}_{\sigma}(t_3)c^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}(t_0)\rangle +\theta(t_1,t_2,t_3,t_0)\langle \hat{c}_{\bar{\sigma}}(t_1)c^{\dagger}_{\bar{\sigma}}(t_2)c^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\rangle \\
-\theta(t_1,t_2,t_0,t_3)\langle \hat{c}_{\bar{\sigma}}(t_1)c^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle +\theta(t_1,t_3,t_0,t_2)\langle \hat{c}_{\bar{\sigma}}(t_1)c^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)c^{\dagger}_{\bar{\sigma}}(t_2)\rangle \\
-\theta(t_2,t_0,t_1,t_3)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle +\theta(t_2,t_1,t_0,t_3)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle \\
-\theta(t_2,t_1,t_3,t_0)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)c^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\rangle +\theta(t_2,t_3,t_1,t_0)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)c^{\dagger}_{\sigma}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\rangle \\
-\theta(t_2,t_3,t_0,t_1)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)c^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\rangle +\theta(t_2,t_0,t_3,t_1)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle \\
-\theta(t_3,t_1,t_0,t_2)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle +\theta(t_3,t_0,t_1,t_2)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle \\
-\theta(t_3,t_2,t_1,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_3),\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\rangle +\theta(t_3,t_2,t_0,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\rangle \\
-\theta(t_3,t_0,t_2,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\rangle +\theta(t_3,t_1,t_2,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\rangle
\end{split}
\end{align}
\end{widetext}
In what follows, we group terms by the order of the operators $\hat{c},\hat{c}^{\dagger}$ to help with the book-keeping. There are four terms with each of the six orderings, \begin{enumerate}[label=(\roman*)]
    \item $\hat{c}\hat{c}\hat{c}^{\dagger}\hat{c}^{\dagger}$
    \item $\hat{c}\hat{c}^{\dagger}\hat{c}^{\dagger}\hat{c}$
    \item $\hat{c}\hat{c}^{\dagger}\hat{c}\hat{c}^{\dagger}$
    \item $\hat{c}^{\dagger}\hat{c}\hat{c}\hat{c}^{\dagger}$
    \item $\hat{c}^{\dagger}\hat{c}\hat{c}^{\dagger}\hat{c}$
    \item $\hat{c}^{\dagger}\hat{c}^{\dagger}\hat{c}\hat{c}$.
\end{enumerate} In evaluating these expectation values we must take more care, because the operator $\hat{n}_{\bar{\sigma}}$, which appears in the time-dependence of $\hat{c}_{\sigma}(t)$ and $\hat{c}_{\sigma}^{\dagger}(t)$, no longer commutes with everything when we have operators which depend on $\bar{\sigma}$. In each of the terms, only one of the four Fermionic states $\ket{0},\ket{\sigma},\ket{\bar{\sigma}},\ket{\sigma\bar{\sigma}}$ gives a nonzero contribution to the trace. The terms are as follows:
\begin{widetext}
\begin{enumerate}[label=(\roman*)]
    \item $\hat{c}\hat{c}\hat{c}^{\dagger}\hat{c}^{\dagger}$ terms: \\
    \begin{itemize}[label={}]
        \item 1.\underline{$t_0>t_1>t_2>t_3$}:  \begin{align}
            -\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle =-Ae^{P(-t_0,-t_1,t_2,t_3)} \exp\big[i\big(F(t_1)-F(t_2)\big)\big]
        \end{align}
        \item 2. \underline{$t_1>t_0>t_2>t_3$}:  \begin{align}
            \langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle =Ae^{P(-t_1,-t_0,t_2,t_3)} \exp\big[i\big(-F(t_2)+F(t_0)\big)\big]
        \end{align}
        \item 3. \underline{$t_1>t_0>t_3>t_2$}: \begin{align}
            -\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle = -Ae^{P(-t_1,-t_0,t_3,t_2)} \exp[i\big(F(t_0)-F(t_3)\big)]
        \end{align}
        \item 4.\underline{$t_0>t_1>t_3>t_2$}: \begin{align}
            \langle \hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle =Ae^{P(-t_0,-t_1,t_3,t_2)}\exp[i\big(F(t_1)-F(t_3)\big)]
        \end{align}
    \end{itemize}
   \item $\hat{c}\hat{c}^{\dagger}\hat{c}^{\dagger}\hat{c}$:
   \begin{itemize}[label={}]
       \item 5.\underline{$t_0>t_2>t_3>t_1$}: \begin{align}\begin{split}
           &-\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)=-Ae^{P(-t_0,t_2,t_3,-t_1)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]\exp\big[i\big(-F(t_2)+F(t_0)\big)\big]\end{split}
       \end{align}
       \item 6.\underline{$t_0>t_3>t_2>t_1$}: \begin{align}\begin{split}
           \langle &\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\rangle =Ae^{P(-t_0,t_3,t_2,-t_1)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]\exp\big[i\big(-F(t_3)+F(t_0)\big)\big]
           \end{split}
       \end{align}
       \item 7.\underline{$t_1>t_3>t_2>t_0$}: \begin{align}\begin{split}
           &-\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\rangle =- Ae^{P(-t_1,t_3,t_2,-t_0)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]
           \exp\big[i\big(-F(t_3)+F(t_1)\big)\big]\end{split}
       \end{align}
       \item 8. \underline{$t_1>t_2>t_3>t_0$}: \begin{align}\begin{split}
      &\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\rangle= Ae^{P(-t_1,t_3,t_2,-t_0)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]\exp\big[i\big(-F(t_2)+F(t_1)\big)\big]\end{split}
       \end{align}
   \end{itemize}
   \item $\hat{c}\hat{c}^{\dagger}\hat{c}\hat{c}^{\dagger}$ terms: \\

   \begin{itemize}[label={}]
       \item 9.\underline{$t_0>t_2>t_1>t_3$}: \begin{align}\begin{split}
           \langle &\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle = Ae^{P(-t_0,t_2,-t_1,t_3)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]\exp\big[i\big(F(t_0)-F(t_2)+F(t_1)-F(t_3)\big)\big]\end{split}
       \end{align}
       \item 10. \underline{$t_0>t_3>t_1>t_2$}: \begin{align}
           -\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle =-Ae^{P(-t_0,t_3,-t_1,t_2)}
       \end{align}
       \item 11. \underline{$t_1>t_2>t_0>t_3$}: \begin{align}
           -\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle =-Ae^{P(-t_1,t_2,-t_0,t_3)}
       \end{align}
       \item 12. \underline{$t_1>t_3>t_0>t_2$}: \begin{align}\begin{split}
           &\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle = Ae^{P(-t_1,t_3,-t_0,t_2)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big] \exp\big[i\big(F(t_1)-F(t_3)+F(t_0)-F(t_2)\big)\big]\end{split}
       \end{align}
   \end{itemize}
   \item $\hat{c}^{\dagger}\hat{c}\hat{c}\hat{c}^{\dagger}$ terms: \\
   \begin{itemize}[label={}]
       \item 13.\underline{$t_2>t_0>t_1>t_3$}: \begin{align}\begin{split}
           &-\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle =-Ae^{P(t_2,-t_0,-t_1,t_3)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]\exp[i\big(F(t_1)-F(t_3)\big)]\end{split}
       \end{align}
       \item 14. \underline{$t_2>t_1>t_0>t_3$}: \begin{align}\begin{split}
          &\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle= Ae^{P(t_2,-t_1,-t_0,t_3)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]\exp[i\big(F(t_0)-F(t_3)\big)]\end{split}
       \end{align}
       \item 15. \underline{$t_3>t_1>t_0>t_2$}: \begin{align}\begin{split}
           &-\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle =-Ae^{P(t_3,-t_1,-t_0,t_2)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big] \exp[i\big(F(t_0)-F(t_2)\big)]\end{split}
       \end{align}
       \item 16.\underline{$t_3>t_0>t_1>t_2$} \begin{align}\begin{split}
           &\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle = Ae^{P(t_3,-t_0,-t_1,t_2)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\\&\times\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]\exp[i\big(F(t_1)-F(t_2)\big)]\end{split}
       \end{align}
   \end{itemize}
   \item $\hat{c}^{\dagger}\hat{c}\hat{c}^{\dagger}\hat{c}$ terms. \\
 \begin{itemize}[label={}]
       \item 17. \underline{$t_2>t_1>t_3>t_0$}: \begin{align}\begin{split}
           &-\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\rangle =-Ae^{P(t_3,-t_0,t_2,-t_1)}\exp\big[\frac{2g^2(t_{min})\beta}{m\omega^2}\big]\exp\big[-\beta U(t_{min})\big]\\&\times\exp\big[2(\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\})\big] \exp[i\big(-F(t_2)+F(t_1)-F(t_3)+F(t_0)\big)]\end{split}
       \end{align}
       \item 18. \underline{$t_2>t_0>t_3>t_1$}: \begin{align}\begin{split}
           &\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle =Ae^{P(t_2,-t_0,t_3,-t_1)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]\end{split}
       \end{align}
       \item 19. \underline{$t_3>t_0>t_2>t_1$}: \begin{align}\begin{split}
           &-\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\rangle =-Ae^{P(t_3,-t_0,t_2,-t_1)}\exp\big[\frac{2g^2(t_{min})\beta}{m\omega^2}\big]\exp\big[-\beta U(t_{min})\big]\\& \times\exp\big[2(\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\})\big]\exp[i\big(-F(t_3)+F(t_0)-F(t_2)+F(t_1)\big)]\end{split}
       \end{align}
       \item 20. \underline{$t_3>t_1>t_2>t_0$}: \begin{align}\begin{split}
           \langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\rangle =Ae^{P(t_3,-t_1,t_2,-t_0)}\exp\big[\frac{g^2(t_{min})\beta}{2m\omega^2}\big]\exp\big[\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\}\big]\end{split}
       \end{align}
   \end{itemize}
   \item $\hat{c}^{\dagger}\hat{c}^{\dagger}\hat{c}\hat{c}$ terms: \\
  \begin{itemize}[label={}]
       \item 21. \underline{$t_2>t_3>t_0>t_1$}: \begin{align}\begin{split}
       &-\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\rangle =-Ae^{P(t_2,t_3,-t_0,-t_1)}\exp\big[\frac{2g^2(t_{min})\beta}{m\omega^2}\big]\exp\big[-\beta U(t_{min})\big]\\& \times\exp\big[2(\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\})\big]\exp[i\big(F(t_1)-F(t_2)\big)]\end{split}
       \end{align}
       \item 22. \underline{$t_2>t_3>t_1>t_0$}: \begin{align}\begin{split}
          &\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_1)\rangle =Ae^{P(t_2,t_3,-t_1,-t_0)}\exp\big[\frac{2g^2(t_{min})\beta}{m\omega^2}\big]\exp\big[-\beta U(t_{min})\big]\\&\times\exp\big[2(\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\})\big] \exp[i\big(F(t_1)-F(t_2)\big)]\end{split}
       \end{align}
       \item 23. \underline{$t_3>t_2>t_1>t_0$}: \begin{align}\begin{split}
           &-\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\rangle =-Ae^{P(t_3,t_2,-t_1,-t_0)}\exp\big[-\beta U(t_{min})\big]\exp\big[\frac{2g^2(t_{min})\beta}{m\omega^2}\big]\\&\times\exp\big[2(\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\})\big] \exp[i\big(F(t_0)-F(t_3)\big)]\end{split}
       \end{align}
       \item 24. \underline{$t_3>t_2>t_0>t_1$}:
       \begin{align}\begin{split}
           &\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\rangle =Ae^{P(t_3,t_2,-t_0,-t_1)}\exp\big[\frac{2g^2(t_{min})\beta}{m\omega^2}\big]\exp\big[-\beta U(t_{min})\big]\\&\times\exp\big[2(\beta\mu-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)+C(t_3)-C(t_1)\})\big]\exp\big[i\big(F(t_1)-F(t_3)\big)]\end{split}
       \end{align}
   \end{itemize}
\end{enumerate}
\end{widetext}
This fully determines the two-particle Green's function in the mixed spin case. This, combined with the expression found for when the spins are the same, determines the two-particle time ordered Green's function of the local Holstein-Hubbard model.

\section{Introduction} \label{sec:intro}
Pump-probe spectroscopy is one of the most important tools for studying complex materials with strongly correlated degrees of freedom. The femtosecond time-resolution allows one to examine excitation and relaxation phenomena from different degrees of freedom (electronic, spin, lattice) at their natural timescales. The technique has been used to probe many different materials in the time-domain: the classification of charge-density-wave insulators \cite{CDWinsulators,CDWPumpProbe} and charge-stripe materials \cite{epcOscillations}; examining the charge gap of spin-density-wave materials \cite{SDWPumpProbe}; determining the electron-phonon coupling in a metal \cite{EPCinAmetal}; studying the response in high $T_c$ cuprates \cite{Tc1, Tc2,Tc3,PDHSSC} and in iron-based superconductors \cite{iron}. In a recent experimental study \cite{iron}, evidence of bosonic mode coupling was seen in LiFeAs. The high $T_c$ cuprates are also relevant systems, since the degree of electron-phonon coupling involved in superconductivity in these materials is not completely understood. For example, a recent study \cite{2020eph} indicates the superconducting behavior cannot be understood in terms of electron-electron or electron-phonon coupling alone, but only by considering both effects.

Another class of materials that this theory might be relevant for are those that have been observed to have polarons (or possibly bipolarons) generated by nonequilibrium pumps. It is difficult to experimentally differentiate between polarons and bipolarons, so some of these systems may be creating bipolarons, which naturally occur in the atomic limit and often are strongly localized. In Ref.~\onlinecite{Polaron1}, the authors use pump-probe spectroscopy to measure polaron lifetimes in a mixed-valence perovskite material under high pressures. They also argue that photoexcitation is analogous to applying a high pressure to this system, evidenced by the response of the phonons and the fact that they both lead to an intervalence charge transfer. More recently, experiments show a photoinduced a transition to a long-lived polaronic state in lead halide perovskite \cite{Polaron2}, and in Ref.~\onlinecite{Polaron3} the authors demonstrate a technique for qualitatively estimating the magnitude and shape of a photoinduced polaronic distortion via femtosecond diffuse x-ray scattering techniques.

Improving the capabilities of pump-probe experiments is currently an active field of research \cite{2019techniqueProbe,3ShotPP2020, PPSpec}, and as time and energy resolution increase, we can expect to obtain a clearer picture of the behavior of electron-phonon coupling, electron-electron coupling, and the interplay between them, in complex materials and in nonequilibrium situations.

The pump is an ultrafast and high-power electric field pulse, which may selectively excite the electronic or lattice degrees of freedom in a material. The subsequent probe is often a second electric field pulse, which reveals information about the material after it has been excited by the initial pulse. By varying the time delay between the pump and probe pulses, experiments can ``watch'' as a material relaxes back to equilibrium. A common feature of these experiments is that the transient states exhibit properties not found in equilibrium.

Another interesting question is whether the pump pulse dynamically modifies either the electron-phonon coupling or the electron-electron coupling in pump-probe experiments. The Shen group  demonstrated\cite{EPCEvidence} that the electron-phonon coupling in the high $T_c$ cuprates is related to a kink in the electronic dispersion near the Fermi energy. More recent work concludes that the strength of the kink is determined by a combination of electron-phonon and electron-electron coupling \cite{huang2020KINK}. A softening of the kink has been observed in pump-probe experments\cite{EPCquenching}, but such softening alone is not sufficient to establish a dynamic change in the electron-phonon coupling due to nonequilibrium effects \cite{JimLexKink}. A recent experimental study of electron thermalization in laser-excited nickel\cite{RecentNickelPaper}, suggests a reduction of electron repulsion may explain the observation of an increasing spectral red shift with increasing electron temperature.

The ability to directly extract the electron-phonon coupling from experiment is an area of active interest that has seen exciting recent developments. Namely, a new analysis technique called frequency-domain angle-resolved photo-emission spectroscopy (FD-ARPES) enables observation of the electron-phonon coupling directly from experimental data \cite{DirectEPC}, and has been used to measure the electron-phonon coupling during a photo-induced insulator to metal transition \cite{DirectEPCexperiment}.

In this work, we allow for arbitrary changes in both the electron-electron and electron-phonon coupling without determining their time-dependence from a self-consistent reaction to the pump pulse. Similar studies have been done previously \cite{Werner_2013,Werner_2015} for the lattice model, but we are not aware of work that varies both electron-electron and electron-phonon couplings as functions of time. We also note that the single and two particle atomic Green's functions presented here are the required inputs to a systematic strong-coupling-expansion on the lattice, generalizing previous work \cite{JimSCE}, which applied only to the Hubbard model.
Recently,\cite{SCE2020} a strong-coupling-expansion of the Holstein-Hubbard model with constant coupling parameters found a rich landscape of phases, including distinct pair-density-waves, spin-density-waves, charge-density-waves, and superconducting regions as well as regions of phase separation.

In this work, we provide an exact (equilibrium and nonequilibrium) solution to the atomic Hubbard-Holstein model for arbitrary electron-electron and electron-phonon couplings. In addition to serving as the starting point for strong-coupling-based perturbation theories, it also provides an important benchmark for nonequilibrium calculations, where exact results are rare. The reason why this problem can be solved exactly lies in the fact that the electron number for each spin is conserved in the atomic limit. But, even with that, the dynamics is sufficiently complex that one needs to use advanced techniques to factorize the evolution operator in order to determine the exact results presented here.

The summary of the paper is as follows. In Sec.~\ref{sec:formalism}, we define the time-dependent Hamiltonian and provide general expressions for the spectral moments, which are analyzed in Sec~\ref{sec:moms}. In Sec.~\ref{sec:timedep} we evaluate the time evolution operator, the partition function, and the time dependence of the field operators. These allow us to calculate the single-particle Green's function in Sec.~\ref{sec:onepart}, verify the moments in Sec.~\ref{sec:moms2} and calculate the two-particle Green's function in Sec.~\ref{sec:twopart}. Our conclusions are presented in Sec.~\ref{sec:conclusion}. Appendices summarize some of the more cumbersome formulas.

\section{Formalism} \label{sec:formalism}
 The Holstein-Hubbard model is the simplest model Hamiltonian which describes both electron-electron and electron-phonon interactions \cite{Hubbard,Holstein}. It includes electrons that hop from one lattice site to another and Einstein phonons with the same frequency on each site. The electrons interact with the phonons via a direct coupling of the charge to the phonon coordinate. The electrons also repel each other via an on-site Coulomb repulsion. In the atomic limit, the electrons do not hop. This means we cannot represent electric-field effects with vector potentials and these systems never have electric current flow. The Holstein Hubbard model for the atomic site then becomes
 \begin{equation}
    \hat{\mathcal H}(t)=\hbar\omega\left (\hat{a}^\dagger\hat{a}+\frac12\right )+\big(g(t)\hat{x}-\mu\big)(\hat{n}_{\uparrow}+\hat{n}_{\downarrow})+U(t)\hat{n}_{\uparrow}\hat{n}_{\downarrow}
\end{equation}
where $\hat{n}_{\sigma}=\hat{c}^{\dagger}_{\sigma}\hat{c}_{\sigma}$ is the fermionic number operator, $\hat{c}_{\sigma}^{\dagger}~(\hat{c}_{\sigma})$ create (destroy) an electron with spin $\sigma$, $\mu$ is the chemical potential, $\hat{a}^{\dagger}~(\hat{a})$ are the bosonic creation (annihilation operators), $\hat{x}=\sqrt{\frac{\hbar}{2m\omega}}(\hat{a}+\hat{a}^{\dagger})$ is the phonon coordinate operator (with $\omega$ the phonon frequency), $g(t)$ is the electron-phonon coupling, and $U(t)$ is the electron-electron coupling. We drop the zero point energy from the Hamiltonian because it just adds an overall constant shift to all energies and we set $\hbar=1$ for the remainder of this work. Note that because $[\hat{\mathcal H}(t),\hat{n}_{\sigma}]=0$, the total electron number (and the electron number for each spin) is conserved. This symmetry is what allows us to derive an exact expression for the retarded Green's function; it is immediately broken once hopping between lattice sites is introduced. Note that we start our system in equilibrium in the limit as $t\to-\infty$ and all of the time dependence in the coupling constants starts at some finite time after $t_{min}$.

The retarded Green's function is defined as
\begin{equation}
\label{gfdef}
    g^R_{\sigma}(t_1,t_2)=-i\theta(t_1-t_2)\big\langle \big\{\hat{c}_{\sigma}(t_1),\hat{c}_{\sigma}^{\dagger}(t_2)\big\}_+\rangle
\end{equation}
where $\theta(t)$ is
the Heaviside function, $\langle \dots\rangle$ indicates the thermal average at the initial time $\frac{1}{\mathcal Z}\text{Tr}\big\{e^{-\beta\hat{\mathcal H}(t_{min})} \dots\big\}$, and all time-dependence of operators is in the Heisenberg representation,
\begin{equation} \hat{O}_H(t)=\hat{U}^{\dagger}(t,t_{min})\hat{O}\hat{U}(t,t_{min})\end{equation} with $\hat{U}(t,t_{min})$ the time-evolution operator, \begin{equation}
    \hat{U}(t,t_{min}) = \mathcal{T}_t\exp\left[-i\int_{t_{min}}^t\!\!\!dt'\ \hat{\mathcal H}(t')\right].
\end{equation}
Here, $\mathcal{T}_t$ is the time-ordering operator, which places the latest times to the left and the Hamiltonian is in the Schr\" odinger representation; $t_{min}$ is a time where the system is in thermal equilibrium and is before any coupling constants start varying with time. The spectral moments of the Green's function are defined via
\begin{equation}
    \mu^{Rn}_{\sigma}(t_{ave})=-\frac{1}{\pi}\int_{-\infty}^{\infty}\,d\omega\ \omega^n\text{Im}\,g^R_{\sigma}(t_{ave},\omega)
\end{equation}
where
\begin{eqnarray}
    &~&g_{\sigma}^R(t_{ave},\omega)=\nonumber\\
    &~&~~\int_0^{\infty}\,dt_{rel}\ e^{i\omega t_{rel}}g^R_{\sigma}\left (t_{ave}+\frac{t_{rel}}{2},t_{ave}-\frac{t_{rel}}{2}\right )
\end{eqnarray}
and we are now using Wigner coordinates of average and relative time,  defined by $t_{ave}=\frac{t_{1}+t_{2}}{2}$ and $t_{rel}=t_1-t_2$. The moments can be rewritten as derivatives with respect to relative time via
\begin{align}
\label{derivs}
    &\mu_{\sigma}^{Rn}(t_{ave})=\lim_{t_{rel}\to 0}\\
    &~~~\text{Re}\bigg\langle i^{n}\frac{d^n}{dt^n_{rel}}\left\{\hat{c}_{\sigma}\left (t_{ave}+\frac{t_{rel}}{2}\right ),\hat{c}_{\sigma}^{\dagger}\left (t_{ave}-\frac{t_{rel}}{2}\right )\right\}_+\bigg\rangle.\nonumber
\end{align}
 Finally we define the time-ordered two-particle atomic Green's function by
 \begin{equation}
 \label{2partgf}
    \mathcal{G}_{\sigma\sigma'}(t_0,t_1,t_2,t_3)=-\langle \mathcal{T}_t \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma'}(t_1)\hat{c}_{\sigma'}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle.
\end{equation}
We will solve for both of these correlation functions below.

\section{Spectral Moment Sum Rules} \label{sec:moms}
 By repeated application of Heisenberg's equation of motion, the spectral moment sum rules can be written as commutators of the field operators with the Hamiltonian plus terms that arise from derivatives of time-dependent coupling constants in the Hamiltonian. After some extensive algebra, we find that the first four moments satisfy \cite{moment2,moments1,erratum,erratum2}
\begin{equation}
    \mu^{R0}_{\sigma}(t_{ave})=\text{Re}\bigg\langle\big\{\hat{c}_{\sigma}(t_{ave}),\hat{c}_{\sigma}^{\dagger}(t_{ave})\big\}_+\bigg\rangle,
\end{equation}
\begin{equation}
    \mu_{\sigma}^{R1}(t_{ave})=\text{Re}\bigg\langle \big\{\big[\hat{c}_{\sigma}(t_{ave}),\hat{\mathcal H}(t_{ave})\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\big\}_+\bigg\rangle,
\end{equation}
\begin{widetext}
\begin{equation}
    \mu_{\sigma}^{R2}(t_{ave})=\text{Re}\bigg\langle \bigg\{\big[[\hat{c}_{\sigma}(t_{ave}),\hat{\mathcal H}(t_{ave})],\hat{\mathcal H}(t_{ave})\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle,
\end{equation}
and
\begin{align}
\begin{split}
    \mu_{\sigma}^{R3}(t_{ave})=\text{Re}\bigg\langle \bigg\{\bigg[\big[[\hat{c}_{\sigma}(t_{ave}),\hat{\mathcal H}(t_{ave})],\hat{\mathcal H}(t_{ave})\big],\hat{\mathcal H}(t_{ave})\bigg],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle\\ +\text{Re}\frac{i}{2}\bigg(\bigg\langle\bigg\{\big[[\hat{c}_{\sigma}(t_{ave}),\frac{\partial \hat{\mathcal H}(t_{ave})}{\partial t_{ave}}],\hat{\mathcal H}(t_{ave})\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle +\bigg\langle\bigg\{\big[[\hat{c}_{\sigma}(t_{ave}),\hat{\mathcal H}(t_{ave})],\frac{\partial \hat{\mathcal H}(t_{ave})}{\partial t_{ave}}\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle\bigg)\\ -\frac{1}{4}\text{Re}\bigg\langle\big[\hat{c}_{\sigma}(t_{ave}),\frac{\partial^2 \hat{\mathcal H}(t_{ave})}{\partial t_{ave}^2}\big],\hat{c}_{\sigma}^{\dagger}(t_{ave})\bigg\}_+\bigg\rangle.
    \end{split}
\end{align}
\end{widetext}
Evaluating the commutators in the atomic limit, we obtain the following:
\begin{equation}
    \mu_{\sigma}^{R0}(t_{ave})=1,
\end{equation}
\begin{equation}
\label{mom1}
    \mu_{\sigma}^{R1}=-\big(\mu-U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle -g(t_{ave})\langle \hat{x}(t_{ave})\rangle\big),
\end{equation}
\begin{align}
\label{mom2}
    \mu_{\sigma}^{R2}(t_{ave})&=\big(\mu-U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle -g(t_{ave})\langle \hat{x}(t_{ave})\rangle\big)^2\\
    &+U(t_{ave})^2
    \big[\langle \hat{n}_{\bar{\sigma}}\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^2\big]\nonumber\\
    &+g^2(t_{ave})\big[\langle \hat{x}^2(t_{ave})\rangle-\langle \hat{x}(t_{ave})\rangle^2\big]\nonumber\\
    &+2U(t_{ave})g(t_{ave})\big[\langle \hat{n}_{\bar{\sigma}}\hat{x}(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle \langle \hat{x}(t_{ave})\rangle\big],\nonumber
\end{align}
\begin{widetext}
\begin{align}
\label{mom3}
    \mu_{\sigma}^{R3}(t_{ave})&=-\big(\mu-g(t_{ave})\langle \hat{x}(t_{ave})\rangle -U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle\big)^3+g^3(t_{ave})\big(\langle \hat{x}^3(t_{ave})\rangle -\langle \hat{x}(t_{ave})\rangle^3\big)\\&-3\mu g^2(t_{ave})\big(\langle \hat{x}^2(t_{ave})\rangle-\langle \hat{x}(t_{ave})\rangle^2\big)
    +U^3(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^3\big)-3\mu U^2(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^2\big)\nonumber\\
    &-6\mu g(t_{ave})U(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\hat{x}(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle\langle \hat{x}(t_{ave})\rangle\big)+3U^2(t_{ave})g(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\hat{x}(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^2\langle \hat{x}(t_{ave})\rangle\big)\nonumber\\
    &+3U(t_{ave})g^2(t_{ave})\big(\langle \hat{n}_{\bar{\sigma}}\hat{x}^2(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle\langle \hat{x}(t_{ave})\rangle^2 \big)\nonumber\\
    &+g(t_{ave})\omega^2\langle \hat{x}(t_{ave})\rangle +\frac{g^2(t_{ave})}{2m}(2\langle \hat{n}_{\bar{\sigma}}\rangle +1) -\frac{1}{4}\text{Re}\bigg(\frac{d^2g(t_{ave})}{dt_{ave}^2}\langle \hat{x}(t_{ave})\rangle +\frac{d^2U(t_{ave})}{dt_{ave}^2}\langle \hat{n}_{\bar{\sigma}}\rangle \bigg)\nonumber
\end{align}
\end{widetext}
where all operators are in the Heisenberg representation. These expressions will be employed as an independent check on the correctness of the single-particle Green's function that we derive next.

\section{Time Dependence}
\label{sec:timedep}
\subsection{Time Evolution Operator}
To obtain the Green's functions of Eqs.~(\ref{gfdef}) and  (\ref{2partgf}), we must first determine the time evolution operator. To do this, we start by switching to the interaction representation with respect to the phonon part of the Hamiltonian and following the technique of Ref.~\onlinecite{jimHolstein}. In the interaction representation, we factorize the time evolution operator into two pieces via
\begin{equation}
    \hat{U}(t,t_{min})=e^{-i\hat{a}^\dagger\hat{a}\omega (t-t_{min})}\mathcal{T}_t\exp\left [-i\int_{t_{min}}^t\!\!\!dt^\prime\hat{\mathcal H}_I(t^\prime)\right ]
\end{equation}
where the interaction Hamiltonian has two pieces $\hat{\mathcal H}_I=\hat{\mathcal H}_I^{el-ph}+\hat{\mathcal H}_I^{coul}$. The electron-phonon coupled piece acquires an additional time dependence through the time dependence of the phonon raising and lowering operators
\begin{equation}
    \hat{\mathcal H}^{el-ph}_I(t)=\frac{g(t)}{\sqrt{2m\omega}}\big(e^{i\omega t} \hat{a}^{\dagger}+e^{-i\omega t}
    \hat{a}\big)\big(\hat{n}_{\sigma}+\hat{n}_{\bar{\sigma}}\big),
\end{equation}
while the Coulomb piece is unaltered ($\hat{\mathcal H}_I^{coul}(t)=U(t)\hat{n}_\uparrow\hat{n}_\downarrow$). These two pieces commute with each other, so the time evolution can be further factorized into a piece from the Coulomb part of the Hamiltonian and a piece from the electron-phonon part. The Coulomb piece is also simple, because it requires no time ordering, and is given by $\exp\left (-i\hat{n}_\uparrow\hat{n}_\downarrow\int_{t_{min}}^t\!\!\!dt^\prime U(t^\prime)\right  )$.
The electron-phonon piece is more complicated because the electron-phonon interaction Hamiltonian at two different times does not commute with itself:
\begin{equation}
    \left [\hat{\mathcal H}^{el-ph}_I(t),\hat{\mathcal H}^{el-ph}_I(t')\right ]=\frac{g(t)g(t^\prime)}{m\omega}i\sin\omega(t-t^\prime) (\hat{n}_{\uparrow}+\hat{n}_{\downarrow})^2,
\end{equation}
but the commutator does commute with $\hat{\mathcal H}^{el-ph}_I(t^{\prime\prime})$. This allows us to solve for the time-evolution operator following the strategy used by Landau and Lifshitz\cite{landau_lifshitz} and by Gottfried\cite{gottfried} in solving the driven simple harmonic oscillator.  We first form the function
\begin{equation}
    \hat{w}_I(t,t_{min})=\int_{t_{min}}^t\,dt'\ \hat{\mathcal H}^{el-ph}_I(t'),
\end{equation}
which satisfies the following equal-time commutator:
\begin{eqnarray}
    \left [\hat{\mathcal H}^{el-ph}_I(t),\hat{w}_I(t,t_{min})\right ]&=&i\frac{g(t)}{m\omega}\int_{t_{min}}^t\!\!\!dt^\prime g(t^\prime)\sin\omega(t-t^\prime)\nonumber\\
    &\times&(\hat{n}_{\uparrow}+\hat{n}_{\downarrow})^2.
\end{eqnarray}
This commutator commutes with all other operators in the different terms in the Hamiltonian.
When these two conditions hold, we have the following operator identity:
\begin{align}
          &e^{i\hat{w}_I(t,t_{min})}\left(i\frac{\partial}{\partial t}-\hat{\mathcal H}^{el-ph}_{I}(t)\right )e^{-i\hat{w}_I(t,t_{min})} \nonumber\\
    &~~~~~=i\frac{\partial}{\partial t}-\frac{i}{2}\left [\hat{w}_I(t,t_{min}),\hat{\mathcal H}^{el-ph}_{I}(t)\right ].
\end{align}
This allows us to evaluate the following time-ordered product directly as
\begin{align}
    &\mathcal{T}_t\exp\left[-i\int_{t_{min}}^t\!\!\!dt'\ \hat{\mathcal H}^{el-ph}_{I}(t')\right ]=\exp\left [-i\hat{w}_I(t,t_{min})\right ]\nonumber\\
&~~~~~\times\exp\left(\frac{1}{2}\int_{t_{min}}^t\!\!\!dt'\ \left [\hat{w}_I(t',t_{min}),\hat{\mathcal H}^{el-ph}_{I}(t')\right ]\right).
\end{align}
Applying this result to our specific Hamiltonian, yields the full time-evolution operator
\begin{align}
     &\hat{U}(t,t_{min})=\exp\left[-i\int_{t_{min}}^t\!\!\!dt'\left (-\mu \hat{n}_f+\omega \hat{a}^{\dagger}\hat{a}+\hat{n}_{\sig}\hat{n}_{\bar{\sigma}}U(t')\right )\right ]\nonumber\\ &\times\exp\left[-i\int_{t_{min}}^t\!\!\!dt' \frac{g(t')}{\sqrt{2m\omega}}(e^{i\omega t'}\hat{a}^{\dagger}+e^{-i\omega t'}\hat{a})\hat{n}_f\right]\\ &\times\exp\left[\frac{i}{2m\omega}\int_{t_{min}}^t\!\!\!dt' \int_{t_{min}}^{t'}\!\!\!dt''\ g(t')g(t'')\sin\omega(t'-t'')\hat{n}_f^2\right]
     \nonumber
\end{align}
where $\hat{n}_f=\hat{n}_{\uparrow}+\hat{n}_{\downarrow}$ is the total electron number operator. Note that because $\hat{n}_{\sigma}^2=\hat{n}_{\sigma}$, as required by the Pauli exclusion principle, we also have that $\hat{n}_f^2=(\hat{n}_{\uparrow}+\hat{n}_{\downarrow})^2=\hat{n}_{\uparrow}+\hat{n}_{\downarrow}+2\hat{n}_\uparrow \hat{n}_\downarrow$.
\subsection{Partition Function}
Calculating thermal averages, requires us to also determine the atomic partition function, which is defined to be
\begin{align}
    \mathcal{Z}&=\text{Tr}_{b,f}\big\{\hat{U}(t=-i\beta+t_{min},t_{min})\big\}\nonumber\\&=\text{Tr}_{b,f}\big\{e^{-\beta\hat{\mathcal H}(t_{min})}\big\},
\end{align}
where the trace is over product states composed of the tensor product of the Fermionic states $\ket{0},\ket{\uparrow},\ket{\downarrow},\ket{\uparrow\downarrow}$ with the Bosonic simple harmonic oscillator number states $|n\rangle=\frac{1}{\sqrt{n!}}\left (\hat{a}^\dagger\right )^n|0\rangle$. Evaluating the trace by using standard methods, we find
\begin{align}
\label{partition}
    \mathcal{Z}&=\Big(n(\omega)+1\Big)\left\{1+2e^{\beta\mu}\exp\left[\frac{\beta g^2(t_{min})}{2m\omega^2}\right]\right.\nonumber\\
    &+\left .\exp\left[\beta(2\mu-U(t_{min}))\right]\exp\left[\frac{2\beta g^2(t_{min})}{m\omega^2}\right]\right\}
\end{align}
with the Bose-Einstein distribution function given by
\begin{align}
\label{bose}
    n(\omega)=\frac{1}{e^{\beta\omega}-1}.
\end{align}
\subsection{Time Dependence of Field Operators}
Next we need the time dependence of the Fermionic field operators in the Heisenberg representation,
\begin{equation}
    \hat{c}_{\sigma}(t)=\hat{U}^{\dagger}(t,t_{min})\hat{c}_{\sigma}\hat{U}(t,t_{min})
\end{equation}
and similarly for the conjugate \begin{equation}
    \hat{c}_{\sigma}^{\dagger}(t)=\hat{U}^{\dagger}(t,t_{min})\hat{c}_{\sigma}^{\dagger}(t)\hat{U}(t,t_{min}).
\end{equation}
To evaluate these operator expressions, we note that the only operator in $\hat{U}(t,t_{min})$ that does not commute with $\hat{c}_{\sigma}$ is $\hat{n}_{\sigma}$. Hence, we need to compute
\begin{equation}
    f(\alpha)=e^{i\alpha \hat{n}_{\sigma}}\hat{c}_{\sigma}e^{-i\alpha \hat{n}_{\sigma}}.
\end{equation}
Differentiating gives us
\begin{equation}
    \frac{df(\alpha)}{d\alpha}=e^{i\alpha \hat{n}_{\sigma}}[\hat{n}_{\sigma},\hat{c}_{\sigma}]e^{-i\alpha \hat{n}_{\sigma}} =-if(\alpha),
\end{equation}
which then can be integrated to yield
\begin{equation}
    f(\alpha)=e^{-i\alpha}\hat{c}_{\sigma}.
\end{equation}\\
Hence, we find that\\
\begin{widetext}
\begin{align}
\hat{c}_{\sigma}(t)&=\exp\left[-i\int_{t_{min}}^t\!\!\!dt'\left (-\mu +\hat{n}_{\bar{\sigma}} U(t')\right )-\frac{i}{\sqrt{2m\omega}}\int_{t_{min}}^t\!\!\!dt'\ g(t')\left (e^{i\omega t'}\hat{a}^{\dagger}+e^{-i\omega t'}\hat{a}\right )\right]\nonumber\\&\times\exp\left[\frac{i}{2m\omega}\int_{t_{min}}^t\!\!\!dt'\ \int_{t_{min}}^{t'}\!\!\!dt''\ g(t')g(t'')\sin\omega(t'-t'')(1+2\hat{n}_{\bar{\sigma}})\right]\hat{c}_{\sigma}
\end{align}
and for the conjugate
\begin{align}
\hat{c}_{\sigma}^{\dagger}(t)&=\exp\left[i\int_{t_{min}}^t\!\!\!dt'\ \left (-\mu+\hat{n}_{\bar{\sigma}}U(t')\right ) +\frac{i}{\sqrt{2m\omega}}\int_{t_{min}}^t\!\!\!dt'\ g(t')\left (e^{i\omega t'}\hat{a}^{\dagger}+e^{-i\omega t'}\hat{a}\right )\right]\nonumber\\&\times\exp\left[-\frac{i}{2m\omega}\int_{t_{min}}^t\!\!\!dt'\ \int_{t_{min}}^{t'}\!\!\!dt''\ g(t')g(t'')\sin\omega(t'-t'')(1+2\hat{n}_{\bar{\sigma}})\right]\hat{c}_{\sigma}^{\dagger}.
\end{align}
\end{widetext}
With these time-dependent operators we are able to exactly evaluate the single- and two-particle Green's functions.

\section{Single-Particle Green's Function}
\label{sec:onepart}
To evaluate the single-particle Green's function in Eq.~(\ref{gfdef}), we must take a trace over the Hilbert space composed of a direct product of the Fermionic states with the simple harmonic oscillator states.
Newton's generalized binomial theorem enables us to sum the infinite degrees of freedom of the harmonic oscillator to obtain a closed form expression for the Bosonic trace. This again follows the technique used in Ref.~\onlinecite{jimHolstein} to evaluate the partition function. We outline this procedure briefly here, as it is used repeatedly in our derivations.  We have
\begin{align}
    g^R_{\sigma}(t_1,t_2)&=-\frac{i\theta(t_1-t_2)}{\mathcal Z}\nonumber\\
    &\times\text{Tr}_{b,f}\bigg\{e^{-\beta\hat{\mathcal H}(t_{min})}\{\hat{c}_{\sigma}(t_1),\hat{c}_{\sigma}^{\dagger}(t_2)\}_+\bigg\}
\end{align}
The Bosonic trace is calculated in the following way. First, we note that the single Bosonic operators always appear as summed pairs, which can be re-expressed by employing the Weyl form of the Baker-Campbell-Hausdorff identity as
\begin{equation}
e^{\alpha \hat{a}+\beta \hat{a}^{\dagger}}=e^{\alpha \hat{a}}e^{\beta \hat{a}^{\dagger}}e^{-\frac{1}{2}\alpha\beta[\hat{a},\hat{a}^{\dagger}]}.
\end{equation}
Using this expression, we can reduce the Bosonic trace to  constants (with respect to the Bosonic operators) multiplied by the form
\begin{equation}
    \text{Tr}_b\left\{e^{-\beta\omega \hat{a}^{\dagger}\hat{a}}e^{A\hat{a}}e^{B\hat{a}^{\dagger}}\right\},
\end{equation}
where $A$ and $B$ are constants \emph{with respect to the Bosonic operators}. Note that the simple harmonic oscillator states $\ket{n}$ are eigenvectors of the Bosonic number operator $\hat{n}=\hat{a}^{\dagger}\hat{a}$ with eigenvalue $n$, so we have
\begin{align}
     &\text{Tr}_b\left\{e^{-\beta\omega \hat{a}^{\dagger}\hat{a}}e^{A\hat{a}}e^{B\hat{a}^{\dagger}}\right\}=\sum_{n=0}^{\infty}e^{-\beta\omega  n}\langle n|e^{A\hat{a}}e^{B\hat{a}^{\dagger}}|n\rangle.
\end{align}
By expanding each exponential inside the expectation value in a power series, using the properties of the simple harmonic oscillator states and summing the resulting series with Newton's generalized binomial theorem, we evaluate the trace to be
\begin{align}
     &\text{Tr}_b\big\{e^{-\beta\omega \hat{a}^{\dagger}\hat{a}}e^{A\hat{a}}e^{B\hat{a}^{\dagger}}\big\}\nonumber\\
     &~~~~~=\frac{1}{1-e^{-\beta\omega}}\exp\left[\frac{1}{1-e^{-\beta\omega}}AB\right].
\end{align}
Using this result allows us to evaluate the Bosonic trace in the definition of the Green's function, which becomes
\begin{align}
\label{bosetrace}
    &g_{\sigma}^R(t_1,t_2)=\frac{-i\theta(t_1-t_2)\left(n(\omega)+1\right)}{\mathcal Z}e^{i\mu(t_1-t_2)}\nonumber\\ &~~~\times\exp\left[\left(-\frac{1}{2}-n(\omega)\right)\big|C(t_1)-C(t_2)\big|^2\right]\nonumber\\ &~~~\times
    \text{Tr}_f\bigg\{\exp\left[\frac{\beta g^2(t_{min})\hat{n}_f^2}{2m\omega^2}+\beta\mu \hat{n}_f-\beta U(t_{min})\hat{n}_{\uparrow}\hat{n}_{\downarrow}\right]\nonumber\\
    &~~~\times\exp\left[-i\int_{t_2}^{t_1}\,dt\ U(t)\hat{n}_{\bar{\sigma}}\right]\nonumber\\
    &~~~\times\exp\left[i(1+2\hat{n}_{\bar{\sigma}})\big(I(t_1)-I(t_2)\big)]\right]\nonumber\\
   &~~~\times \exp\left[\frac{2ig(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\big\{C(t_1)-C(t_2)\big\}\hat{n}_f\right]\nonumber\\
    &~~~\times\left\{e^{i\text{Im}\{C^{*}(t_1)C(t_2)\}}\hat{c}_{\sigma}\hat{c}_{\sigma}^{\dagger}+ e^{-i\text{Im}\{C^{*}(t_1)C(t_2)\}}\hat{c}_{\sigma}^{\dagger}\hat{c}_{\sigma}\right\}\bigg\}
\end{align}
where
\begin{equation}
\label{Coft}
    C(t)=\frac{1}{\sqrt{2m\omega}}\int_{t_{min}}^t\!\!\!dt'\ g(t') e^{i\omega t'},
\end{equation}
and \begin{align}
\label{Ioft}
    I(t)=\frac{1}{2m\omega}\int_{t_{min}}^t\!\!\!dt' \int_{t_{min}}^{t'}\!\!\!dt''\ g(t')g(t'')\sin\omega(t'-t'')
\end{align}
are functions we define to make this result more readable.
We can now evaluate the Fermionic trace over the four possible states, with the result
\begin{widetext}
\begin{align}
\label{final1}
g_{\sigma}^R(t_1,t_2)&=-\frac{i\theta(t_1-t_2)\big (n(\omega)+1\big)}{\mathcal Z}e^{i\mu(t_1-t_2)}
\exp\bigg[\left(-\frac{1}{2}-n(\omega)\right)\big|C(t_1)-C(t_2)\big|^2\bigg]\nonumber\\&\times
\bigg\{e^{i\text{Im}\{C^{*}(t_1)C(t_2)\}}\bigg(\exp\big[i\big(I(t_1)-I(t_2)\big)\big]\nonumber\\
&+\exp\left[\beta(\mu+\frac{g^2(t_{min})}{2m\omega^2})\right] \exp\left[-i\int_{t_2}^{t_1}\,dt'\ U(t')\right]\exp\big[3i(I(t_1)-I(t_2))\big]
\exp\left[\frac{2ig(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{C(t_1)-C(t_2)\}\right ]\bigg)\nonumber\\
&+e^{-i\text{Im}\{C^{*}(t_1)C(t_2)\}}\bigg( \exp\left[\beta(\frac{\mu+g^2(t_{min})}{2m\omega^2})\right] \exp\big[i\big(I(t_1)-I(t_2)\big)\big] \exp\left[\frac{2ig(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{C(t_1)-C(t_2)\}\right]\nonumber\\
&+\exp\left[\frac{2\beta g^2(t_{min})}{m\omega^2}+\beta(2\mu-U(t_{min}))\right]\exp\left[-i\int_{t_2}^{t_1}\,dt\ U(t)\right]\exp\big[3i(I(t_1)-I(t_2))\big]\nonumber\\
&\times\exp\left[\frac{4ig(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{C(t_1)-C(t_2)\}\right]\bigg)\bigg\}
\end{align}
\end{widetext}

This result is quite complicated. But note that it is an exact expression, so it can be employed to describe the response of Holstein-Hubbard systems with arbitrary time-dependent interactions. It is remarkable that such a complex system has such a compact exact expression for the nonequilibrium Green's functions. We reiterate again that there are two keys to allowing this to occur: first the number operators for each spin are conserved and do not depend on time and second, the commutator of the electron-phonon interaction with itself at two different times is an operator that commutes with all remaining operators in the system. Nevertheless, the final result is not simple. It represents the complex dynamics of electron-phonon coupled systems in the atomic limit. In contrast, the atomic limit Green's functions for the Hubbard model are quite simple.

While it would be nice to have an independent check to verify that these results are correct, we are not aware of any way to do so for the full functions. But we do demonstrate below that this Green's function does satisfy the first four spectral moment sum rules. The first three are evaluated analytically by taking derivatives of Eq.~(\ref{bosetrace}), before we take the Fermionic trace, as in Eq.~(\ref{derivs}). The third moment is verified by taking numerical derivatives of the equation Eq.~(\ref{final1}). We also numerically check the first three moments to ensure we have taken the trace correctly in going from Eq.~(\ref{bosetrace}) to Eq.~(\ref{final1}). While this does not test everything, it is a quite stringent test of accuracy of the final results.

\section{Verification of Moments}
\label{sec:moms2}
\begin{figure}[b]
    \centering
    \includegraphics[width=0.4\textwidth,angle=270]{thirdmomHH.jpg}
    \caption{Comparison of the exact third moment  [given by Eq.~(\ref{mom3})] and the numerical third derivative of Eq.~(\ref{final1}) for non-constant electron-phonon and electron-electron coupling parameters {$g(t)=1-0.05e^{-0.01t}$} and {$U(t)=1-0.01t^2$}.}
    \label{fig:Third}
\end{figure}
Using the definition of the Green's function and Eq.~(\ref{derivs}), allows us to write the moments as derivatives of the causal Green's function with respect to relative time  \begin{equation}
    \mu_{\sigma}^{Rn}(t_{ave})=-\text{Im}\big\{i^n\frac{\partial^n}{\partial t_{rel}^n}g^R_{\sigma}(t_{ave}+\frac{t_{rel}}{2},t_{ave}-\frac{t_{rel}}{2})\big\}\big|_{t_{rel\to 0}}.
\end{equation}
Using this form, we easily find from Eq.~\ref{final1} that $\mu_{\sigma}^{R0}(t_{ave})=1$. By evaluating traces similar to those we evaluated to determine the partition function and the Green's function, we also find exact expressions for the time-dependent expectation values, which appear in the moments of Eqs.~(\ref{mom1}-\ref{mom3}). Following Ref.~\onlinecite{jimHolstein}, we have following two expectation values:
\begin{widetext}
    \begin{equation}
    \langle \hat{x}(t)\rangle =-\langle \hat{n}_f\rangle \left (\frac{g(t_{min})}{m\omega^2}\cos\big(\omega (t-t_{min})\big)+\text{Re}\left\{ie^{-i\omega t}\int_{t_{min}}^t\!\!\!dt'\ \frac{g(t')}{m\omega}e^{i\omega t'}\right\}\right)
\end{equation}
and \begin{equation}
\label{nf}
    \langle \hat{n}_f\rangle =\frac{2e^{\beta\mu}\exp\left[\frac{\beta g^2(t_{min})}{2m\omega^2}\right]+2e^{\beta(2\mu-U(t_{min}))}\exp\left[\frac{2\beta g^2(t_{min})}{m\omega^2}\right]}{1+2e^{\beta\mu}\exp\left[\frac{\beta g^2(t_{min})}{2m\omega^2}\right]
    +\exp\left[\beta(2\mu-U(t_{min}))\right]\exp\left[\frac{2\beta g^2(t_{min})}{m\omega^2}\right]}
\end{equation}
\end{widetext}
The first is time dependent, while the second is time independent (because electron number is conserved in the model).
The remaining expectation values needed to evaluate the analytic expressions for the moments are given in Appendix~\ref{sec:a}. By taking derivatives of Eq.~(\ref{bosetrace}) and simplifying we find \begin{equation}
    \mu_{\sigma}^{R1}=-\mu+U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle +g(t_{ave})\langle \hat{x}(t_{ave})\rangle,
\end{equation}
in agreement with Eq.~(\ref{mom1})) and
\begin{widetext}
\begin{align}
    \mu_{\sigma}^{R2}(t_{ave})&=\big(\mu-U(t_{ave})\langle \hat{n}_{\bar{\sigma}}\rangle -g(t_{ave})\langle \hat{x}(t_{ave})\rangle\big)^2+U(t_{ave})^2
    \big[\langle \hat{n}_{\bar{\sigma}}\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle^2\big]\nonumber\\
    &+g^2(t_{ave})\big[\langle \hat{x}^2(t_{ave})\rangle-\langle \hat{x}(t_{ave})\rangle^2\big]+2U(t_{ave})g(t_{ave})\big[\langle \hat{n}_{\bar{\sigma}}\hat{x}(t_{ave})\rangle-\langle \hat{n}_{\bar{\sigma}}\rangle \langle \hat{x}(t_{ave})\rangle\big],
\end{align}
\end{widetext}
in agreement with Eq.~(\ref{mom2}), proving that the first three moments of the single-particle Green's function hold exactly. For the third moment, taking derivatives by hand becomes too cumbersome, so we resort to taking numerical derivatives to verify that the third moment holds. In Fig.~\ref{fig:Third}, we plot this comparison in the case where $t_{min}=0$, $U(t)=1-0.01t^2$ and $g(t)=1-0.05\exp[-0.1t]$. One can clearly see the third moment is verified numerically. While we plot just one example, this agreement held for all combinations of parameters that we tested, providing strong evidence that the third-moment sum rule is also satisfied by the single-particle Green's function in Eq.~(\ref{final1}).
While this exercise is not definitive, it is hard to imagine how one can satisfy all of these moment sum rules and still have an error (of course, it could still occur with the higher-order terms, but trying to evaluate higher-order moments is also too cumbersome).

\section{Two Particle Green's Function}
\label{sec:twopart}
For the two particle Green's function, there are no known sum rules to compare our results against, so we no longer have an independent check. The definition of the two-particle Green's function is given in Eq.~(\ref{2partgf}). There are two possibilities for the spin $\sigma'$, $\sigma'=\sigma$ or $\sigma'=\bar{\sigma}$.
\subsection{Same Spin Case}
For the case where the spins are the same, we find a number of terms in the time ordering vanish, because $\hat{c}^2_{\sigma}=\hat{c}_{\sigma}^{\dagger 2}=0$ due to the Pauli exclusion principle. In this case, only 8 of the 24 possible time-orderings are non-vanishing. When we expand the time-ordering operator for the same spin Green's function, we obtain
\begin{widetext}
\begin{align}
\label{samespin}
    \mathcal{G}_{\sigma\sigma}(t_0,t_1,t_2,t_3)&=-\langle {\mathcal T}_t\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}(t_1)\hat{c}^{\dagger}_{\sigma}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle
    \\&=\theta(t_0,t_2,t_1,t_3)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle-\theta(t_0,t_3,t_1,t_2)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\rangle \nonumber\\
    &+\theta(t_1,t_3,t_0,t_2)\langle \hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\rangle -\theta(t_1,t_2,t_0,t_3)\langle \hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle \nonumber\\
    &+\theta(t_2,t_0,t_3,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\rangle -\theta(t_2,t_1,t_3,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\rangle \nonumber\\
    &+\theta(t_3,t_1,t_2,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\rangle -\theta(t_3,t_0,t_2,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_2)\hat{c}_{\sigma}(t_1)\rangle. \nonumber
\end{align}

The expectation values are evaluated in the same fashion as we did for the single-particle case. For example, the first term, which is nonzero when $t_0>t_2>t_1>t_3$, is given by
\begin{align}
    &\langle \hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\sigma}(t_2)\hat{c}_{\sigma}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle
    =
    \frac{\big(n(\omega)+1\big)}{\mathcal Z}e^{-i\mu(-t_0+t_2-t_1+t_3)}
    \exp\left[\left(-\frac{1}{2}-n(\omega)\right)\big|C(t_2)-C(t_0)+C(t_1)-C(t_3)\big|^2\right]\\ &~~~~~~\times\exp\big[i\big(I(t_0)-I(t_2)+I(t_1)-I(t_3)\big)\big]
    \exp[P(-t_0,t_2,-t_1,t_3)]\nonumber\\
    &~~~~~~\times\bigg\{1+\exp\left[\beta\left(\mu+\frac{g^2(t_{min})}{2m\omega^2}\right)\right] \exp\left[-i\frac{2g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(t_0)+C(t_2)-C(t_1)+C(t_3)\}\right]\nonumber\\ &~~~~~~\times\exp\left[-i\left(-2I(t_0)+\int_{t_{min}}^{t_0}\!\!\!dt' U(t')+2I(t_2)-\int_{t_{min}}^{t_2}\!\!\!dt'U(t')
    -2I(t_1)+\int_{t_{min}}^{t_1}\!\!\!dt' U(t')+2I(t_3)-\int_{t_{min}}^{t_3}\!\!\!dt' U(t')\right)\right]\bigg\}\nonumber
\end{align}
where $C(t)$ and $I(t)$ are defined in Eqs.~(\ref{Coft}) and (\ref{Ioft}) and $P(t,t', \bar{t},\bar{t}')$ is a prefactor, which results from the repeated application of the Baker-Campbell-Hausdorff theorem, and is defined by the following expression
\begin{align}
        &P((\pm)_1t,(\pm)_2t',(\pm)_3\bar{t},(\pm)_4\bar{t}')=
        \exp\bigg[-i\text{Im}\big\{(\pm)_1(\pm)_2C^{*}(t)C(t')
        +(\pm)_3(\pm)_4C^{*}(\bar{t})C(\bar{t}')\\&~~~~~~+(\pm)_1(\pm)_3C^{*}(t)C(\bar{t})
        +(\pm)_1(\pm)_4C^{*}(t)C(\bar{t}')+(\pm)_2(\pm)_3C^{*}(t')C(\bar{t})
        +(\pm)_2(\pm)_4C^{*}(t')C(\bar{t}')\big\}\bigg]\nonumber
\end{align}
\end{widetext}
where the $(\pm)_i$ is a $+$ if the time $i$ is associated with a $c^{\dagger}$ operator, $-$ for a $c$ operator.
The expressions for other time orderings are given in appendix~\ref{sec:b}.
\subsection{Different Spin Case}
In the case where $\sigma'=\bar{\sigma}$, none of the terms in the time-ordering vanish. Then the Green's function for the mixed spin case is given by
\begin{widetext}
\begin{align}
\label{differentspin}
\mathcal{G}_{\sigma\bar{\sigma}}(t_0,t_1,t_2,t_3)&=
-\theta(t_0,t_1,t_2,t_3)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle+\theta(t_1,t_0,t_2,t_3)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle \nonumber\\
&-\theta(t_1,t_0,t_3,t_2)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle +\theta(t_0,t_1,t_3,t_2)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\rangle \nonumber\\
&-\theta(t_0,t_2,t_3,t_1)\langle \hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle +\theta(t_0,t_2,t_1,t_3)\langle \hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\rangle \nonumber\\
&-\theta(t_0,t_3,t_1,t_2)\langle \hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle +\theta(t_0,t_3,t_2,t_1)\langle \hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle \nonumber\\
&-\theta(t_1,t_3,t_2,t_0)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}(t_0)\rangle +\theta(t_1,t_2,t_3,t_0)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\rangle \nonumber\\
&-\theta(t_1,t_2,t_0,t_3)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle +\theta(t_1,t_3,t_0,t_2)\langle \hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}^{\dagger}_{\bar{\sigma}}(t_2)\rangle \nonumber\\
&-\theta(t_2,t_0,t_1,t_3)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle +\theta(t_2,t_1,t_0,t_3)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\rangle \nonumber\\
&-\theta(t_2,t_1,t_3,t_0)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\rangle +\theta(t_2,t_3,t_1,t_0)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\rangle \nonumber\\
&-\theta(t_2,t_3,t_0,t_1)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}^{\dagger}_{\sigma}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\rangle +\theta(t_2,t_0,t_3,t_1)\langle \hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\rangle \nonumber\\
&-\theta(t_3,t_1,t_0,t_2)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle +\theta(t_3,t_0,t_1,t_2)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\rangle \nonumber\\
&-\theta(t_3,t_2,t_1,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\sigma}(t_0)\rangle +\theta(t_3,t_2,t_0,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}(t_1)\rangle \nonumber\\
&-\theta(t_3,t_0,t_2,t_1)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\sigma}(t_0)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\bar{\sigma}}(t_1)\rangle +\theta(t_3,t_1,t_2,t_0)\langle \hat{c}_{\sigma}^{\dagger}(t_3)\hat{c}_{\bar{\sigma}}(t_1)\hat{c}_{\bar{\sigma}}^{\dagger}(t_2)\hat{c}_{\sigma}(t_0)\rangle.
\end{align}
For example, the term with $t_0>t_2>t_1>t_3$ is given by
\begin{align}
           &\langle c_{\sigma}(t_0)c_{\bar{\sigma}}^{\dagger}(t_2)c_{\bar{\sigma}}(t_1)c_{\sigma}^{\dagger}(t_3)\rangle = \frac{\big(n(\omega)+1\big)}{\mathcal Z}e^{-i\mu(-t_0+t_2-t_1+t_3)}
    \exp\left[\left(-\frac{1}{2}-n(\omega)\right)\big|C(t_2)-C(t_0)+C(t_1)-C(t_3)\big|^2\right]\nonumber\\
    &~~~~~\times\exp\big[i\big(I(t_0)-I(t_2)+I(t_1)-I(t_3)\big)\big]
           e^{P(-t_0,t_2,-t_1,t_3)}
           \exp\left[\frac{\beta g^2(t_{min})}{2m\omega^2}\right]\exp\big[\beta\mu\big]
           \nonumber\\ &~~~~~\times\exp\left[-2i\frac{g(t_{min})}{\sqrt{2m\omega^3}}\text{Re}\{-C(-t_0)+C(t_2)-C(t_1)+C(t_3)\}\right]\nonumber\\
      &~~~~~\times\exp\bigg[-i\big(-2I(t_0)+\int_{t_{min}}^{t_0}\!\!\!dt' U(t')+2I(t_2)-\int_{t_{min}}^{t_2}\!\!\!dt' U(t')
    -2I(t_1)+\int_{t_{min}}^{t_1}\!\!\!dt' U(t')+2I(t_3)-\int_{t_{min}}^{t_3}\!\!\!dt' U(t')\big)\bigg]\bigg\}
       \end{align}
      \end{widetext}
While there are many terms, the evaluation is again straightforward, but lengthy. Results for the other time orderings in the mixed spin case are also presented in Appendix~\ref{sec:b}.

\section{Conclusion}
\label{sec:conclusion}
In this work, we have found an exact result for the single-particle and two-particle correlation functions of the nonequilibrium Holstein-Hubbard model (with time dependent interactions) in the atomic limit. This includes all of the complex dynamical effects associated with the coupling of electrons to quantum phonons and the possible change in time of the electron-electron and electron-phonon interactions. Such dynamics are, of course, quite complex, so the results are cumbersome. While there is no independent check on their correctness, we did partially verify them by comparing them to exact results of the first four nonequilibrium spectral moments of the single-particle retarded Green's function, finding exact agreement.  We also calculated the two-particle Green's functions for this model. Both of these Green's functions are needed to calculate the lattice Green's functions within a strong-coupling-expansion framework. We believe such future work will eventually  provide insight into the behavior of complex materials in pump-probe spectroscopy experiments, especially those that have bipolarons. In the future, it may be interesting to also expand these results to the case with a time-dependent  phonon frequency and to consider non-linear couplings of the electrons to the phonons. This would be even more impactful if one can properly determine how the change in electron occupancy during a pump modifies the interaction parameters self-consistently. Perhaps these exact results could be employed within a more complicated formalism to do just that for systems that are near the atomic limit. But, even as standalone results, these expressions provide potential benchmarks to nonequilibrium many-body calculations that examine the behavior of strongly coupled electron-phonon systems near the atomic limit.

redpajama_set_name: RedPajamaArXiv