\section{Introduction}

The physics of strongly interacting quantum many-body systems
is, despite decades of research, far from being fully understood. 
Apart from a few exactly solvable models, even their ground state properties 
are a subject of ongoing discussions, see, e.g., \cite{L04,YJNIHY03}.
Even more challenging is the investigation of non-equilibrium 
properties of those interacting many-body systems.

In the limit of infinite dimensions, the Dynamical Mean-field Theory (DMFT) 
has been successfully applied to various problems such as the Mott-Hubbard 
transition \cite{RZK92,RKZ94,GK92,GKKR96,AGPTW10,AGPTW11} or the non-equilibrium dynamics in Mott insulators \cite{EKW09,ATE14,EKW10,WTE12,CWE14,FT06} 
by mapping the system to an effective single site (i.e., zero-dimensional) 
problem. 

Furthermore, one-dimensional strongly interacting systems have been studied 
by employing exact diagonalization, see, e.g., \cite{K16,LHS88,KLA07,KTEG08,RB04}, 
time-dependent Density Matrix Renomalization Group variational techniques (t-DMRG),  
see, e.g., \cite{DKSV04,S11,WF04}, or Jordan-Wigner transformations \cite{JW28,LP75,BO01,BKOSV04}. 
Moreover, exact analytical solutions were used to study non-equilibrium 
dynamics \cite{D05} and ground state properties.
However, the findings for one-dimensional systems cannot be 
easily transferred to higher dimensions.
For example, thermalization in higher dimensions is very different 
from thermalization in one dimension since for the latter
the energy re-distribution cannot occur via two-body collisions
(due to energy and momenta conservation).

For strongly interacting systems in higher dimensions, 
the methods described above run into difficulties.
For example, the generalization of t-DMRG to higher dimensions 
(such as tensor networks) is limited by the exponential scaling 
with the system width \cite{CV09}, especially if correlations spread 
across the system. 
For a very similar reason (exponential scaling of Hilbert space), 
the method of exact diagonalization is limited to 
small system sizes \cite{K16}.
On the other hand, the DMFT is only well understood in the limit of infinite 
dimensions and, due to the mapping to an effective single-site problem, does 
not capture energy and momentum transfer of quasi-particles or long-range correlations 
\cite{WEFWBH18}.

In order to bridge this gap,
we established over the last years a perturbative hierarchical method which 
is valid for large coordination numbers $Z\gg1$.
This method allows for a systematic study of non-equilibrium properties in 
strongly interacting systems in large (but finite) dimensions \cite{NS10,QNS12,QKNS14,KNQS14,NQS14,NQS16,QS19,QS19a}.

The hierarchical expansion is based on a controlled expansion 
of the $n$-point reduced density matrices into correlated parts.
At zeroth order in our expansion, we have the single-site density matrix,
\begin{equation}
\hat{\rho}_\mu=\mathrm{Tr}_{\not{\mu}}\{\hat\rho\}=
\mathcal{O}(Z^0)\,,
\end{equation}
where $\mathrm{Tr}_{\not{\mu}}$ denotes the trace over all lattice sites 
but $\mu$. 
The correlated part of the two-site density matrix, 
\begin{equation}
\hat{\rho}^\mathrm{corr}_{\mu\nu}
=\mathrm{Tr}_{\not{\mu}\not{\nu}}\{\hat\rho\}-\hat{\rho}_\mu \hat{\rho}_\nu
\,, 
\end{equation}
is of order $\mathcal{O}(1/Z)$, the three-point correlator 
\begin{equation}
\hat{\rho}^\mathrm{corr}_{\mu\nu\lambda}
=\mathrm{Tr}_{\not{\mu}\not{\nu}\not{\lambda}}
\{\hat\rho\}-\hat{\rho}^\mathrm{corr}_{\mu\nu} 
\hat{\rho}_\lambda-\hat{\rho}_{\mu\lambda}^\mathrm{corr} 
\hat{\rho}_\nu -\hat{\rho}_{\nu\lambda}^\mathrm{corr} 
\hat{\rho}_\lambda-\hat{\rho}_\mu \hat{\rho}_\nu \hat{\rho}_\lambda 
\end{equation}
is of order $\mathcal{O}(1/Z^2)$, and so on.
The time-evolution of the $n$-point correlations is derived 
directly from the von Neumann equations.
With this expansion, we were able to study quenches 
across phase boundaries \cite{NS10}, ground state properties and non-equilibrium 
dynamics of quantum correlations in the bosonic and fermionic 
Hubbard model \cite{QKNS14,KNQS14}.
Although the perturbative scheme is set up for $Z\gg1$, we 
found even qualitative agreement with exact diagonalization  in one dimension \cite{KNQS14}. 

The derivation of the hierarchy was based on the real-time 
evolution of the system's density matrix.
In the following, we shall extend our hierarchical 
approach to double-time correlation functions.
As a concrete example, we will apply the method to the Bose-Hubbard model 
in the Mott insulating phase.
In section \ref{H} we derive the equations of motion 
up to order $\mathcal{O}(1/Z)$. 
A general proof of the hierarchy is given in the Appendix.
In section \ref{BHequil}, we show how Green functions and finite 
temperature correlation functions are related to the hierarichal equations.
In section \ref{BHnonequil1} we study the non-equilibrium dynamics at finite 
temperatures and in section \ref{BHnonequil2} we derive the double-time correlation 
functions for a quantum quench of the Bose-Hubbard model.



\section{The Double-Time Hierarchy}\label{H}

The lattice system under consideration is described by the Bose-Hubbard 
Hamiltonian ($\hbar=1$)
\begin{align}\label{bosehubbard}
\hat{H}=-\frac{J(t)}{Z}\sum_{\mu,\nu}T_{\mu\nu}\hat{b}^\dagger_\mu\hat{b}_\nu+\sum_\mu \left[\frac{U}{2}\hat{n}_\mu(\hat{n}_\mu-1)-\mu_0 
\hat{n}_\mu\right]\,, 
\end{align}
where $\hat b_{\mu}^\dagger$ and $\hat b_{\nu}$ are the bosonic creation 
and annihilation operators at the lattice sites $\mu$ and $\nu$, respectively.
The time-dependent hopping rate is denoted by $J(t)$, where we have 
factored out the coordination number $Z$.
The second term describes the on-site repulsion $U$ with the particle 
number operator $\hat n_\mu$ and the chemical potential~$\mu_0$.


The system (\ref{bosehubbard}) can only be solved exactly for a small number of lattice sites. 
For large systems one relies on suitable approximation schemes.
As mentioned in the introduction, we shall employ 
a hierarchy of correlations for large coordination numbers $Z$.
For arbitrary operators $\hat A_\mu(t)$ and $\hat B_\mu(t)$ 
we have (see the Appendix)
\begin{align}
\langle \hat A_\mu(t) \hat B_\mu(t')\rangle&=\mathcal{O}(Z^0)\\
\langle \hat A_\mu(t) \hat B_\nu(t')\rangle^\mathrm{corr}&=\langle \hat A_\mu(t) \hat B_\nu(t')\rangle-\langle \hat A_\mu(t)\rangle\langle \hat B_\nu(t')\rangle\nonumber\\
&=\mathcal{O}(1/Z),\, \mathrm{etc.}
\end{align}
We assume the system to be in the Mott insulating phase
(with integer filling)
where the on-site repulsion $U$ dominates over the hopping rate $J$. 
To order $\mathcal{O}(Z^0)$, the density matrix of the system 
is given by
\begin{align}
\hat{\rho}=\bigotimes_\mu \hat\rho_\mu, \quad
\hat\rho_\mu=\sum_n p_n(t)|n\rangle_\mu \langle n| 
\end{align}
with the on-site probabilities $p_n(t)$.
We define the operators $\hat{P}^{n,m}_\mu=|n\rangle_\mu\langle m|$
and find for the on-site expectation values
\begin{align}\label{singlesite}
&i\partial_t\langle \hat{P}^{n,n+1}_\mu(t) \hat{P}^{m+1,m}_\mu(t')\rangle\nonumber\\
&=(U n -\mu_0)\langle \hat{P}^{n,n+1}_\mu(t) \hat{P}^{m+1,m}_\mu(t')\rangle
+\mathcal{O}(1/Z)\,.
\end{align}
The two-site correlation functions with $\mu\neq\nu$
evolve according to
\begin{align}\label{twosite}
&\left(i\partial_t-U n +\mu_0\right)\langle \hat{P}^{n,n+1}_\mu(t) 
\hat{P}^{m+1,m}_\nu(t')\rangle^\mathrm{corr}\nonumber\\
=&-\frac{J(t)}{Z}\sum_{\kappa\neq \mu,\nu}T_{\mu\kappa}\sqrt{n+1}[p_n(t)-p_{n+1}(t)]\nonumber\\
&\times \sum_{l=0}^\infty\sqrt{l+1}\langle \hat{P}^{l,l+1}_\kappa(t)\hat{P}^{m+1,m}_\nu(t')\rangle^\mathrm{corr}\nonumber\\
&-\frac{J(t)}{Z}T_{\mu\nu}\sqrt{n+1}[p_n(t)-p_{n+1}(t)]\nonumber\\
&\times \sum_{l=0}^\infty\sqrt{l+1}\langle \hat{P}^{l,l+1}_\nu(t)\hat{P}^{m+1,m}_\nu(t')\rangle+\mathcal{O}(1/Z^2)\,.
\end{align}
In equation (\ref{twosite}) we employed the hierarchy 
in order to separate
three-point expectation values, 
\begin{align}
&\langle\hat{P}_\mu^{n,n}(t) \hat{P}^{l,l+1}_\kappa(t)
\hat{P}^{m+1,m}_\nu(t')\rangle=\nonumber\\
&p_n(t)\langle\hat{P}^{l,l+1}_\kappa(t)
\hat{P}^{m+1,m}_\nu(t')\rangle^\mathrm{corr}+\mathcal{O}(1/Z^2)\,, \mathrm{etc.} 
\end{align}
In the following we shall study only the 
evolution up to order $\mathcal{O}(1/Z)$.
Thus we limit our considerations to the linear evolution 
of quasi-particle excitations.
All results in this work will be an immediate 
consequence of the evolution equations (\ref{singlesite}) and (\ref{twosite}).

Note that the hierarchical expansion can be extended to 
arbitrary high orders in $1/Z$, e.g., the interaction 
among the quasi-particle excitation can be taken into account \cite{QS19a}.
The proof of whole $1/Z$-hierarchy for double-time 
correlation functions involving arbitrary $n$-point 
correlations is given in the appendix \ref{PH}.








\section{Double-time expectation values in equilibrium}\label{BHequil}





The hierarchical set of equations for two-time expectation values 
induces a hierarchy for double-time Green's functions.
In the following we will use them for the computation 
of equilibrium correlation functions in the Bose-Hubbard model at finite temperatures.
A thermal expectation value of an operator $\hat O $ is defined as
\begin{eqnarray}
\langle \hat O\rangle_\mathrm{therm}=\frac{\mathrm{tr}(\hat{O}e^{-\beta\hat{H}})}{\mathrm{tr}
(e^{-\beta\hat{H}})}\,,
\end{eqnarray}
where $\beta=1/k_b T$ is the inverse temperature.
The retarded and advanced Green's functions for 
operators $\hat{A}_\mu(t)$ and $\hat{B}_\nu(t')$  are given by \cite{Z60}
\begin{eqnarray}\label{green}
G^r_{A_\mu,B_\nu}(t,t')&=&\langle\langle \hat{A}_\mu(t);\hat{B}_\nu(t')\rangle\rangle^{r}\nonumber\\
&=&-i \Theta(t-t')
\langle[\hat{A}_\mu(t),\hat{B}_\nu(t')]\rangle
\end{eqnarray}
and 
\begin{eqnarray}
G^a_{A_\mu,B_\nu}(t,t')&=&\langle\langle \hat{A}_\mu(t);\hat{B}_\nu(t')\rangle\rangle^{a}\nonumber\\
& =&i \Theta(t'-t)
\langle[\hat{A}_\mu(t),\hat{B}_\nu(t')]\rangle\,.
\end{eqnarray}
Both Green functions, although they are subject to different boundary 
conditions, obey the same differential equation.
We choose $\hat{A}_\mu=\hat{P}^{n,n+1}_\mu$ and $\hat{B}_\mu=\hat{P}^{m+1,m}_\mu$
and find from equations (\ref{singlesite}) and (\ref{green}) the relation
\begin{align}
&i\partial_t\langle\langle \hat{P}^{n,n+1}_\mu(t); \hat{P}^{m+1,m}_\mu(t')\rangle\rangle
 =\delta(t)\delta_{m,n}(p_n-p_{n+1})\nonumber\\
&+(Un-\mu_0)\langle\langle \hat{P}^{n,n+1}_\mu(t); \hat{P}^{m+1,m}_\mu(t')\rangle\rangle+\mathcal{O}(1/Z)\,.
\end{align}
The on-site probabilities $p_n$ are time-independent since we evaluate 
all quantities in a thermal equilibrium state
and the hopping rate $J$ is assumed to be time-independent.
An analogous equation for $\mu\neq\nu$ follows from equation (\ref{twosite}). 
Via a Fourier transform w.r.t. space and time, we can turn
the differential equations for the Green functions into  a set of algebraic equations,
\begin{align}
&(\omega-U n+\mu_0 )\langle\langle \hat{P}^{n,n+1}; \hat{P}^{m+1,m}\rangle\rangle_\omega=
\frac{\delta_{m,n}}{2\pi}(p_n-p_{n-1})\,,\label{onsitegreen}\\
&(\omega-U n+\mu_0 )\langle\langle \hat{P}^{n,n+1}; \hat{P}^{m+1,m}\rangle\rangle_{\mathbf{k},\omega}=\nonumber\\
&-JT_\mathbf{k}\sqrt{n+1}(p_n-p_{n+1})\sum_{l=0}^\infty \sqrt{l+1}\nonumber\\
&\times\left[\langle\langle \hat{P}^{l,l+1}; \hat{P}^{m+1,m}\rangle\rangle_{\mathbf{k},\omega}-\langle\langle \hat{P}^{l,l+1}; \hat{P}^{m+1,m}\rangle\rangle_{\omega}\right]\label{twositegreen}\,.
\end{align}



Using the spectral decomposition of operators, it can be shown that 
thermal expectation-value of double-time correlation functions and 
the Green's functions are related by \cite{Z60}
\begin{align}\label{zub}
\langle 
\hat{B}(t)\hat{A}(t')\rangle_\mathrm{therm}=&i\lim_{\epsilon\rightarrow0^+}\int_{-\infty}^\infty 
d\omega\frac{e^{i \omega (t-t')}}{e^{\beta \omega}-1}\nonumber \\
\times&\left[\langle\langle \hat{A};\hat{B}\rangle\rangle_{\omega+i\epsilon}-\langle\langle 
\hat{A};\hat{B}\rangle\rangle_{\omega-i\epsilon} \right]\,.
\end{align}
The exact relation (\ref{zub}) permits the calculation of time-independent 
expectation values ($t=t'$) or double-time expectation values 
($t\neq t'$) at finite temperatures.
From equation (\ref{onsitegreen}) we obtain the probabilities for $n$ bosons on a lattice site,
\begin{align}
p_n&=\frac{e^{-\beta\left(\frac{U}{2}n(n-1)-\mu_0 n\right)}}{\sum_{m=0}^\infty 
e^{-\beta\left(\frac{U}{2}m(m-1)-\mu_0 m\right)}}\,.\label{thermalpn}
\end{align}




From equation (\ref{twositegreen}) we find after some algebra the two-site correlator in order $1/Z$,
\begin{align}\label{corrbb}
&\langle \hat{b}^\dagger_\mu(t)\hat{b}_\nu(t')\rangle_\mathrm{therm}=
\frac{\mathrm{Tr}\left(e^{-\beta\hat{H}}\hat{b}^\dagger_\mu(t)\hat{b}_\nu(t')\right)}
{\mathrm{Tr}\left(e^{-\beta\hat{H}}\right)}\nonumber\\
=&\frac{1}{ \pi N}\int d\omega\, \mathcal{P}\frac{1}{e^{\beta \omega}-1}\sum_\mathbf{k}\mathrm{Im}\left(\frac{JT_\mathbf{k}G(\omega+i\epsilon)^2}
{1+JT_\mathbf{k}G(\omega+i\epsilon)}\right)\nonumber\\
&\times e^{i\mathbf{k}\cdot (\mathbf{x}_\mu-\mathbf{x}_\nu)+i\omega(t-t')}\,,
\end{align}
where we introduced the thermal Green function \cite{NQS16}
\begin{align}
G(\omega)&=\sum_{n=0}^\infty \frac{(n+1)(p_n-p_{n+1})}{\omega-Un+\mu_0}\,.
\end{align}
For a hyper-cubic lattice in $D$ dimensions, the 
correlation function (\ref{corrbb}) can be expressed with
modified Bessel functions $I_\alpha(z)$,
\begin{align}
&\langle \hat{b}^\dagger_\mu(t)\hat{b}_\nu(t')\rangle_\mathrm{therm}
=-\frac{1}{\pi}\int d\omega\,\mathcal{P}\frac{1}{e^{\beta \omega}-1} \int_0^\infty ds e^{-s+i\omega(t-t')}\nonumber\\
&\times\mathrm{Im}\left[ G(\omega+i\epsilon) \prod_{j=1}^D 
(-1)^{\triangle \mathbf{x}_{\mu\nu}^j}I_{\triangle \mathbf{x}_{\mu\nu}^j}
\left(\frac{J G(\omega+i\epsilon)s}{D}\right)\right]\,,
\end{align}
where $\triangle \mathbf{x}_{\mu\nu}^j=x_{\mu}^j-x_{\nu}^j$
is the $j$-component of the distance between two lattice sites.
$\triangle \mathbf{x}_{\mu\nu}^j$ is always an integer since we rescaled the distances with the lattice constant. 




In Fig.~\ref{diagcorrdoubletime} we evaluate the double-time 
correlation function for different values of $\beta U$.
The density plots show a temperature-dependent light-cone structure.
For low temperatures it is viable to neglect occupation numbers $n>2$.
The dynamics of the correlation functions w.r.t.
to the time difference $t_-=(t-t')/2$ at unit filling
is determined by the frequencies 
\begin{align}\label{partholetherm}
\Omega_\mathbf{k}=&U-J(1-3p_0) T_\mathbf{k}\nonumber\\
&-\sqrt{U^2-6 J(1-3p_0) U T_\mathbf{k}+J^2(1-3p_0)^2T_\mathbf{k}^2}
\end{align}
with $p_0=p_2=1/(2+e^{\beta U/2})$ being the probability to find zero or two particles on a lattice site.
The maximum group velocity, given by 
$v^\mathrm{max}=\mathrm{max}_\mathbf{k}|\nabla_\mathbf{k} \Omega_\mathbf{k}|$,
determines the light-cone structure of the correlations.
At zero temperature, the propagation velocity of the 
correlations is maximal, see Fig.~\ref{diagcorrdoubletime} (d).
In this limit, only the doublon-holon excitations $\langle P^{10}_\mu(t)P^{12}_\nu(t')\rangle$,  
$\langle P^{21}_\mu(t)P^{01}_\nu(t')\rangle$,  
$\langle P^{10}_\mu(t)P^{01}_\nu(t')\rangle$,
and
$\langle P^{21}_\mu(t)P^{12}_\nu(t')\rangle$
are relevant. 

For finite temperatures, we have $1-3p_0<1$.
According to equation (\ref{partholetherm}), this implies  
a shrinking propagation velocity, see Fig.~\ref{diagcorrdoubletime} (c).
If the temperature is increased even further, the 
approximation of thermal particle-hole excitations
is not valid anymore since also occupation numbers $n>2$ play a role.
Therefore, the analytical result (\ref{partholetherm}) 
does not determine the light-cone structure anymore, 
see Fig.~\ref{diagcorrdoubletime} (a)~and~(b).

The population of higher excited states is also reflected 
in the time-evolution of the correlations at a fixed distance $|x_\mu-x_\nu|$.
At $\beta U=\infty$, only eigen-modes of order $\mathcal{O}(U)$ play a role whereas 
for large temperatures the eigen-modes of order $\mathcal{O}(nU)$,~$n>1$
dominate the time-evolution, see Fig.~\ref{ThermalCorrl1}.

\begin{figure}[h]
\includegraphics[width=6.8cm]{ThermalCorrl1.pdf}
\caption{Double-time correlation for $|x_\mu-x_\nu|=\sqrt{2}$ in two dimensions 
at different temperatures.
The position of the maximum is shifting with temperature 
approximately as $(t_-J)\sim 1/(1-3 p_0)$.
}\label{ThermalCorrl1}
\end{figure}


 \newpage

 
\begin{widetext}

\begin{figure}
\includegraphics[width=17cm]{diagcorrdoubletime.pdf}
\caption{Time-dependence of the two-time correlation functions $|\langle \hat{b}^\dagger(t)\hat{b}(t')\rangle_\mathrm{therm}|^2$ 
along the diagonal for a lattice in two dimensions for unit filling and $J/U=0.1$.
The parameters are (a) $\beta U=1$, (b) $\beta U=3$, (c) $\beta U=5$, (d) $\beta U=\infty$.
The maximum velocity of the correlations decreases with the temperature 
of the thermal ensemble.
When the temperature grows, also highly energetic modes contribute to the dynamics.
The amplitude of the correlations decreases with increasing temperature.
}\label{diagcorrdoubletime}
\end{figure}

\end{widetext}



\section{Equilibration after a quantum quench}\label{BHnonequil1}



In this section we turn to the non-equilibrium dynamics of excitations 
in the Bose-Hubbard system after a quantum quench at finite 
temperatures from $J=0$
to $J>0$.
The evolution of time-local quantities
can be deduced from double-time hierarchy for $t=t'$ 
of from the single-time hierarchy in \cite{NS10}.
Afterwards we compare the pre-thermalized expectation 
values at $t\rightarrow \infty$ with the corresponding thermal expectation values.


We assume the system to be initially in a 
thermal state which is determined by the 
density matrix
\begin{align}
\hat{\rho}^\mathrm{in}=\bigotimes_\mu \hat\rho_\mu, \quad
\hat\rho_\mu=\sum_n p_n|n\rangle_\mu \langle n|\,, 
\end{align}
where the $p_n$ are given by the thermal distribution (\ref{thermalpn}).
The hopping rate is switched suddenly from $J=0$
to a finite value which is still in the 
Mott regime $J<J_\mathrm{crit}$.
The dynamics of the on-site probabilities follow the equations
\begin{align}\label{onsite}
&i\partial_t \langle \hat{P}^{n,n}_\mu\rangle=-\frac{J}{Z}\sum_{\kappa}T_{\mu\kappa}
\sum_m \sqrt{m+1}\\
\times&\bigg[\sqrt{n}\big(\langle \hat{P}_\mu^{n,n-1}\hat{P}_\kappa^{m,m+1}\rangle^\mathrm{corr}
-\langle \hat{P}_\kappa^{m+1,m}\hat{P}_\mu^{n-1,n}\rangle^\mathrm{corr}\big)\nonumber\\
&-\sqrt{n+1}\big(\langle \hat{P}_\mu^{n+1,n}\hat{P}_\kappa^{m,m+1}\rangle^\mathrm{corr}
-\langle \hat{P}_\kappa^{m+1,m}\hat{P}_\mu^{n,n+1}\rangle^\mathrm{corr}\big)\bigg]\nonumber\,,
\end{align}
and the two-point correlations evolve according to
\begin{align}\label{twosites}
&[i\partial_t+U(n-m)]\langle \hat{P}^{n+1,n}_\mu\hat{P}^{m,m+1}_\nu\rangle^\mathrm{corr}\nonumber\\
=-&\frac{J}{Z}\sqrt{n+1}(\langle \hat{P}^{n+1,n+1}_\mu\rangle-\langle \hat{P}^{n,n}_\mu\rangle)\nonumber\\
&\times\sum_{\kappa\neq\mu,\nu}T_{\kappa\mu}
\langle \hat{b}^\dagger_\kappa\hat{P}^{m,m+1}_\nu\rangle^\mathrm{corr}\nonumber\\
-&\frac{J}{Z}\sqrt{m+1}(\langle \hat{P}^{m,m}_\nu\rangle-\langle \hat{P}^{m+1,m+1}_\nu\rangle)\nonumber\\
&\times\sum_{\kappa\neq\mu,\nu}T_{\kappa\nu}
\langle \hat{P}^{n+1,n}_\mu \hat{b}_\kappa\rangle^\mathrm{corr}\nonumber\\
-&\frac{J}{Z}T_{\mu\nu}\sqrt{n+1}\sqrt{m+1}\langle \hat{P}^{n+1,n+1}_\mu\rangle\langle \hat{P}^{m,m}_\nu\rangle\nonumber\\
+&\frac{J}{Z}T_{\mu\nu}\sqrt{n+1}\sqrt{m+1}\langle \hat{P}^{n,n}_\mu\rangle\langle \hat{P}^{m+1,m+1}_\nu\rangle\,.
\end{align}
In the limit of zero temperature, the doublon-holon correlation 
functions decouple from correlation functions which involve 
occupation numbers $n>2$.
This case has been discussed in references \cite{NS10,QKNS14}.
However, we shall see that a finite temperature alters
the pre-thermalization dynamics.



The equations (\ref{onsite}) and (\ref{twosites})
are nonlinear due to the back-reaction of the 
on-site quantities onto the correlations.
This back-reaction can be neglected since 
it is of order $1/Z^2$ and we have 
$\langle \hat P^{n,n}_\mu\rangle=P_n=p_n+\mathcal{O}(1/Z)$.
Exploiting the translational symmetry, 
we can simplify these equations by a spatial 
Fourier transformation
\begin{align}\label{Fourier}
\langle \hat{P}_\mu^{n,n-1} \hat{P}_\nu^{m-1,m}\rangle=\frac{1}{N}\sum_\mathbf{k}f^{n,m}_\mathbf{k}e^{i\mathbf{k}\cdot(\mathbf{x}_\mu-\mathbf{x}_\nu)}\,,
\end{align}
where $N$ denotes the number of lattice sites.
Now, the on-site probabilities $P_n=\langle \hat{P}^{n,n}_\mu\rangle $ evolve according to
\begin{align}\label{pn}
i\partial_t  P_n=&-\frac{J}{N}\sum_{\mathbf{k}}T_\mathbf{k}
\sum_m \sqrt{m+1}\nonumber\\
&\times\bigg[\sqrt{n}\big(f_\mathbf{k}^{n,m+1}
-f_\mathbf{k}^{m+1,n}\big)\nonumber\\
&-\sqrt{n+1}\big(f_\mathbf{k}^{n+1,m+1}
-f_\mathbf{k}^{m+1,n+1}\big)\bigg]\,,
\end{align}
whereas for the Fourier components of the correlations 
we find
\begin{align}\label{bb}
&[i\partial_t+U(n-m)]f_\mathbf{k}^{n+1,m+1}\nonumber\\
&=-JT_\mathbf{k}\sqrt{n+1}(p_{n+1}-p_n)\sum_{s=0}^\infty 
\sqrt{s+1}f_\mathbf{k}^{s+1,m+1}\nonumber\\
&-JT_\mathbf{k}\sqrt{m+1}(p_{m}-p_{m+1})\sum_{s=0}^\infty 
\sqrt{s+1}f_\mathbf{k}^{n+1,s+1}\nonumber\\
&-JT_\mathbf{k}\sqrt{m+1}\sqrt{n+1}(p_{n+1}p_m-p_{n}p_{m+1})\,.
\end{align}
We solved the coupled system of equations 
in two and three dimensions for various temperatures, see Fig.~\ref{dynamics2D3D}.
The on-site probabilities are oscillating around their initial 
thermal value $p_n$.
Especially for high temperatures, the back action of these 
quantities onto the correlation functions can indeed be 
safely neglected since\linebreak $|P_n(t)-P_n(0)|/P_n(0)\ll1$.


In order to gain some analytical insight 
into the equilibration process, we
neglect $p_n$ for $n>2$.
As a consequence, we obtain the particle-hole 
symmetry $p_0=p_2$.
When the system is quenched from a thermal state, 
the oscillations decrease with 
time and the on-site probabilities 
approach for $t\rightarrow\infty $
\begin{align}\label{pnequil}
&P_{0/2,\mathrm{equil}}= p_{0}+\frac{1}{N}\sum_\mathbf{k}
\frac{4J^2T_\mathbf{k}^2(1-4p_0+3p_0^2)}{\omega_\mathbf{k}^2}
\end{align}
with
\begin{align}
\omega_\mathbf{k}=&
\sqrt{\big[U^2-6JT_\mathbf{k}(1-3p_0)U+J^2T_\mathbf{k}^2(1-3p_0)^2\big]}\,.
\end{align}
For the lattice-site correlations we for $t\rightarrow\infty $ the asymptotic 
expression
\begin{align}\label{correquil}
&\langle \hat{b}^\dagger_\mu\hat{b}_\nu\rangle_\mathrm{equil}=
\frac{1}{N}\sum_\mathbf{k}e^{i\mathbf{k}\cdot(\mathbf{x}_\mu-\mathbf{x}_\nu)}
\frac{4JT_\mathbf{k}U(1-4p_0+3p_0^2)}{\omega_\mathbf{k}^2}\,.
\end{align}
Given equations (\ref{pnequil}) and (\ref{correquil}), we conclude that the quench-induced 
change of these quantities becomes smaller with increasing 
temperature.

From the eigen-modes
we can also estimate the maximum propagation 
speed of the expanding correlations.
As before, we can estimate the maximum propagation 
velocity from $v^\mathrm{max}=\mathrm{max}_\mathbf{k}|\nabla_\mathbf{k}\Omega_\mathbf{k}|$.
In a hyper-cubic lattice in $D$ dimensions 
with small $J/U$ we have  
$v^\mathrm{max}=J(3-9p_0)/D$ along the lattice axes
and $v^\mathrm{max}=J(3-9p_0)/\sqrt{D}$ along the diagonals.
As before, the maximum propagation speed becomes smaller 
for increasing temperatures, see Fig.~\ref{vmax}.



 \begin{center}
\begin{figure}[h]
\includegraphics[width=7.3cm]{vmax}
\caption{The maximum velocity of the correlations after a quantum quench 
from $J=0$ to $J/U=0.1$ depends on the temperature of the initial 
ensemble.}\label{vmax}
\end{figure}
\end{center}



\begin{widetext}
 
\begin{figure}
\includegraphics[width=18cm]{dynamics2D3D}
\caption{Equilibration dynamics in two and three dimensions.
The amplitude of the oscillations decreases with increasing temperature 
of the initial state.
For all initial temperatues, the system prethermalizes faster in three dimensions than it does 
in two dimensions.
In our computation, we did not restrict the occupation number.
Nevertheless, we see that also for finite but sufficiently low temperatures, the 
particle-hole symmetry is approximately valid.
}\label{dynamics2D3D}
\end{figure}

\end{widetext}



The pre-thermalized state is a quasi-stationary state which is not 
necessarily thermal whereas real thermalization occurs on much longer 
time-scales \cite{BBS03,BBW04,KWE11}.
However, we can compare the pre-thermalized state 
with a corresponding thermal state that has the 
same temperature $T$ and filling $\langle \hat{n}_\mu\rangle $.
In order to do so, we compute the asymptotic values 
of the on-site probabilities $P_{n,\mathrm{equil}} $ and the correlation functions 
$\langle \hat{b}_\mu^\dagger\hat{b}_\nu\rangle_\mathrm{equil}$
from equations (\ref{pn}) and (\ref{bb}) without further approximations.
The thermal correlation $\langle \hat{b}_\mu^\dagger\hat{b}_\nu\rangle_\mathrm{therm}$ 
can be deduced from equation (\ref{corrbb}).
For the order $\mathcal{O}(1/Z)$-contribution of the $P_{n,\mathrm{therm}} $ we take
the expression from reference \cite{NQS16}.
In this article, we evaluated the partition function 
of the Bose-Hubbard system up to $\mathcal{O}(1/Z)$ and derived 
various thermal quantities.
Explicitly we have (see equation 21 in \cite{NQS16})
\begin{align}
P_{n,\mathrm{therm}}=&p_n+\frac{1}{N}\sum_\mathbf{k}\sum_{l=-\infty}^{\infty}
\frac{J T_\mathbf{k}}{1+JT_\mathbf{k}G(\omega_l)}\bigg[p_n G(\omega_l)\nonumber\\
&-\frac{(n+1)p_n}{\omega_l -U n+\mu}+\frac{np_n}{\omega_l-U(n-1)+\mu}\nonumber\\
&+\frac{n(p_{n-1}-p_n)}{\beta(\omega_l-U(n-1)+\mu)^2)^2}\nonumber\\
&-\frac{(n+1)(p_n-p_{n+1})}{\beta(\omega_l-U n+\mu)^2}\bigg]
\end{align}
with the bosonic Matsubara modes $\omega_l=2\pi i l/\beta$.


At zero temperature, the values 
differ roughly by a factor of 2 which has been 
stated elsewhere \cite{MK09}.
For sufficiently low temperatures, 
the particle-hole symmetry is valid, $P_0\sim P_2$, which 
can be deduced from comparing Fig.~\ref{thermalasym2D3D} (a)
and (b) or Fig.~\ref{thermalasym2D3D} (d) and (e).
The energy which is transferred 
by a quench to the system becomes 
more and more irrelevant for increasing 
temperatures.
Thus, the difference between thermal and pre-thermalized 
quantities reduces with growing $T$.
From Fig.~\ref{thermalasym2D3D} (c) and (f) we see 
how the increasing temperature diminishes the correlations
between lattice sites. 
However, the thermal next neighbor correlation 
function has its global maximum at finite $T>0$.


 \newpage


\begin{widetext}
 


 

\begin{figure}
\includegraphics[width=17cm]{thermalasym2D3D}
\caption{After a quantum quench at a temperature $T $, the occupation probabilities and the correlation functions approach a prethermalized value (red curves).
This is compared with the corresponding thermal expressions (black curves).
For low temperatures, the onsite-quantities differ roughly by a factor of 2. 
The difference between prethermalized and thermal 
values diminishes with increasing temperatures.
Note that deviations from the particle-hole symmetry at high temperatures are larger
in three dimensions (see (d) and (e)) than in two dimensions (see (d) and (e)).
In (c) and (f) we depict the temperature-dependence of the prethermalized and thermal expectation values of the next neighbour correlations $\langle \hat{b}_\mu^\dagger\hat{b}_\nu\rangle$.
Note that the thermal correlation function has its global maximum at $T>0$.
}\label{thermalasym2D3D}
 \end{figure}


\end{widetext}
\section{Quantum quench and Double-time correlations}\label{BHnonequil2}



A quantum quench changes also the light-cone structure of the double-time correlation functions.
We assume that before the quantum quench the system is in an equilibrium 
state with a finite hopping rate $J_-$.
At time $t=0$, the hopping rate is quenched to a value 
$J_+<J_\mathrm{crit}$.
For simplicity, we restrict our considerations 
to zero temperature.
Analogous to equation (\ref{Fourier}) we define the 
Fourier transform of the two-time correlations
\begin{align}
\langle \hat{P}_\mu^{n,n-1}(t) \hat{P}_\nu^{m-1,m}(t')\rangle=\frac{1}{N}\sum_\mathbf{k}f^{n,m}_\mathbf{k}(t,t')e^{i\mathbf{k}\cdot(\mathbf{x}_\mu-\mathbf{x}_\nu)}\,.
\end{align}
The dynamics of the correlation functions follows directly from equation (\ref{twosite})
and a corresponding relation which takes into account the time-evolution w.r.t.~variable $t'$.
Neglecting the back-reaction onto the on-site probabilities, the Fourier coefficients
evolve up to order $1/Z$ according to
\begin{align}\label{odet1}
i\partial_t f^{11}_\mathbf{k}&= -J^\pm T_\mathbf{k}(f_\mathbf{k}^{11}+\sqrt{2}f_\mathbf{k}^{21})-J^\pm T_\mathbf{k}\,,\\
i\partial_t f^{12}_\mathbf{k}&= -J^\pm T_\mathbf{k}(f_\mathbf{k}^{12}+\sqrt{2}f_\mathbf{k}^{22})\,,\\
i\partial_t f^{21}_\mathbf{k}&= \sqrt{2}J^\pm T_\mathbf{k}(f_\mathbf{k}^{11}+\sqrt{2}f_\mathbf{k}^{21})-U f_\mathbf{k}^{21}+\sqrt{2}J^\pm T_\mathbf{k}\,,\\
i\partial_t f^{22}_\mathbf{k}&= \sqrt{2}J^\pm T_\mathbf{k}(f_\mathbf{k}^{12}+\sqrt{2}f_\mathbf{k}^{22})-U f_\mathbf{k}^{22}\,,
\end{align}
and
\begin{align}\label{odet2}
i\partial_{t'} f^{11}_\mathbf{k}&= J^\pm T_\mathbf{k}(f_\mathbf{k}^{11}+\sqrt{2}f_\mathbf{k}^{12})+J^\pm T_\mathbf{k}\,,\\
i\partial_{t'} f^{12}_\mathbf{k}&= -\sqrt{2}J^\pm T_\mathbf{k}(f_\mathbf{k}^{11}+\sqrt{2}f_\mathbf{k}^{12})+U f_\mathbf{k}^{12}-\sqrt{2}J^\pm T_\mathbf{k}\,,\\
i\partial_{t'} f^{21}_\mathbf{k}&= J^\pm T_\mathbf{k}(f_\mathbf{k}^{21}+\sqrt{2}f_\mathbf{k}^{22})\,,\\
i\partial_{t'} f^{22}_\mathbf{k}&= -\sqrt{2}J^\pm T_\mathbf{k}(f_\mathbf{k}^{21}+\sqrt{2}f_\mathbf{k}^{22})+U f_\mathbf{k}^{22}
\,.
\end{align}
For $t,t'<0$, the correlation functions can be immediately obtained from the ground state 
correlations at $J=J_-$ \cite{QKNS14}.
Due to the quantum quench, the correlation function is not anymore homogeneous in 
time but depends on the relative time $t_-=(t-t')/2$
and on the central time $t_+=(t+t')/2$.
For $t_+ \pm t_-<0$ and $t_+ \pm t_->0$ we find correlations before and after the hopping quench, respectively. 
Correlations between an creation (annihilation) of a 
particle before the quench and an annihilation (creation)
of particle can be obtained for $t_+ + t_-<0$ and 
$t_+ - t_->0$ ($t_+ + t_->0$ and 
$t_+ - t_-<0$). 


The light-cone structure is primarily determined by three different velocities.
For $t,t'<0$, the maximum velocity $v^\mathrm{max}_{--}=\mathrm{max}_\mathbf{k}|\nabla_\mathbf{k} \Omega^{--}_\mathbf{k}|$ can be estimated from 
\begin{align}
\Omega_\mathbf{k}^{--}=U-J_-T_\mathbf{k}-\sqrt{U^2-6 J_- UT_\mathbf{k}+(J_-T_\mathbf{k})^2}\,. 
\end{align}
For $t>0$ and $t'<0$ or $t<0 $ and $ t'>0$, the spread is determined by
both hopping rates $J_-$ and $J_+$.
Thus we have the eigen-modes
\begin{align}
\Omega_\mathbf{k}^{+-}=&\frac{1}{2}\left(U-J_-T_\mathbf{k}-\sqrt{U^2-6 J_- UT_\mathbf{k}+(J_-T_\mathbf{k})^2}\right)\nonumber\\
+&\frac{1}{2}\left(U-J_+T_\mathbf{k}-\sqrt{U^2-6 J_+ UT_\mathbf{k}+(J_+T_\mathbf{k})^2}\right)\,.
\end{align}
from which one can obtain $v^\mathrm{max}_{+-}=\mathrm{max}_\mathbf{k}|\nabla_\mathbf{k} \Omega^{+-}_\mathbf{k}|$.
Finally, if $t,t'>0$, the spread of correlations is dominantly determined by $J_+$.
Here the maximum velocity $v^\mathrm{max}_{++}=\mathrm{max}_\mathbf{k}|\nabla_\mathbf{k} \Omega^{++}_\mathbf{k}|$
can be derived from
\begin{align}
\Omega_\mathbf{k}^{++}=U-J_+T_\mathbf{k}-\sqrt{U^2-6 J_+ UT_\mathbf{k}+(J_+T_\mathbf{k})^2}\,. 
\end{align}
The kink in the light-cone structure due to the quantum quench 
is illustrated in Fig.~\ref{doubletimecorr}.
The complete light-cone structure for the double-time correlators 
in arbitrary dimensions is evaluated in appendix \ref{doubtime}.

\begin{figure}
\includegraphics[width=7cm]{doubletime}
\caption{Double-time correlation function $\langle \hat{b}^\dagger_\mu (t)\hat{b}_\nu (t')\rangle_{T=0}$
in one dimension for a quench at $t=0$ from $J_-/U=0.05$ to $J_+/U=0.15$.
(a): $t_+U=100$, the lightcone has a kink at $t_-=t_+$ where the maximum velocity 
changes from $v^\mathrm{max}_{++}\approx0.40$ to $v^\mathrm{max}_{+-}\approx0.24$.
(b): $t_+U=-100$, the lightcone has a kink at $t_-=-t_+$ where the maximum velocity 
changes from $v^\mathrm{max}_{--}\approx 0.10$ to $v^\mathrm{max}_{+-}\approx0.24$.
}\label{doubletimecorr}
\end{figure}


\section{Conclusions}

We studied equilibrium properties and non-equilibrium 
dynamics of the Bose-Hubbard model in the Mott insulating phase.
To this end, we extended the hierarchy for large coordination 
numbers $Z$ presented in \cite{NS10} to a hierarchy for double-time 
correlation functions.

This enabled us to derive thermal correlation functions 
and we studied the spread of two-time correlation functions at 
finite initial temperatures, see Figs.~\ref{diagcorrdoubletime} and \ref{ThermalCorrl1}.
We restricted our considerations to first order in $1/Z$ which 
governs the free quasi-particle evolution.
As usual, the effective light-cone structure was obtained via a
saddle-point approximation.

As we demonstrated above, the dynamics of strongly interacting 
quantum many-body systems and their thermal properties can be 
both accessed within our approach.
The phenomenon of pre-thermalization after a quantum quench 
was generalized from $T=0$ \cite{QKNS14} to finite temperatures, 
see Figs.~\ref{dynamics2D3D} and \ref{vmax}.
Note that the dynamical equations in order $\mathcal{O}(1/Z)$ 
cover only the short-time evolution of the excitations.
In order to include the long-time evolution which 
is determined by quasi-particle scattering, 
also higher order correlations have to be taken into account.
For real thermalization processes, for example, one has to 
evaluate the hierarichal equations up to order 
$\mathcal{O}(1/Z^3)$ in order to obtain a Boltzmann equation 
\cite{QS19a}.

Apart from these theoretical investigations, our approach can 
also be applied in order to interpret experimental results.  
Since the double-time hierarchy induces a Green function 
hierarchy, the latter can used for the interpretation of 
pump-probe experiments \cite{}.
An excitation which is created by a pump beam at time $t$ 
and measured by a probe beam at time $t'$ is naturally 
described in terms of the lesser Green functions \cite{ATE14,FMSSD17,FKP17}.
Moreover, the Green function hierarchy allows for a controlled
perturbative calculation of the self-energy or the 
susceptibility of strongly correlated systems.

As an outlook, the study of out-of-time-correlators (OTOC), 
a common measure for quantum chaos, should be feasible within our approach \cite{LO69,MSS16}.
In this context one has to extend the double-time hierarchy of correlations to
a hierarchy of multi-time correlations.
This should be possible in complete analogy to the approach 
presented in this work and will be the subject of further studies. 


\acknowledgments


This work was funded by DFG, grant \# 278162697 (SFB 1242) and 398912239.



