Let's assume we have a model, which is initially defined by the tight-binding Hamiltonian with a random on-site energy $f_n$, as follows:
$$H^i=-J\sum_n^{L-1}\left(a_n^\dagger a_{n+1}+h.c\right)+\sum_n^Lf_na_n^\dagger a_n$$
The above single-particle Hamiltonian can be diagonalized in real space with basis change:
$$\eta_m^\dagger=\sum_n U_{nm}a_m^\dagger$$
So, we have $H^i=\sum_m\epsilon_m\eta_m^\dagger\eta_m$, and the many-body ground state can be constructed by creating a particle up to a filling point:
$$\mid\Psi_0^i \rangle=\Pi_m^N\eta_m^\dagger\mid0\rangle$$
where $\mid 0 \rangle$ is the vaccum state, and $N\le L$ is the number of particles.
Now, we do a global quench and turn the random on-site term off:
$$H^f=-J\sum_n^{L-1}\left(a_n^\dagger a_{n+1}+h.c\right)$$
This Hamiltonian is also diagonalizable in both momentum and real space, so all energy and corresponding eigenstates are available. From now on, the dynamic of the model is driven by $H^f$. We are interested in calculating the unequal time correlation defined by:
$$\left<\Psi_o^i|a_r(t)a_{r'}^\dagger(0)|\Psi_o^i\right>=\left<\Psi_o^i|e^{iH^ft}a_r(0)e^{-iH^ft}a_{r'}^\dagger(0)|\Psi_o^i\right>$$
The question is, how to calculate numerically the above equation?
Best Answer
For fermionic Gaussian states, we can use the covariance matrix to completely characterize the state, $$ C_{mn} = \langle \Psi|a^{\dagger}_{m} a_{n}|\Psi\rangle. $$ This also applies to mixed state. The correlation matrix is a semi-positive Hermitian matrix, with eigenvalues between 0 and 1.
Given a set of occupying orbitals labelled by $$ \eta^{\dagger}_k = \sum_{m=1}^L U_{km}a^{\dagger}_m, $$ where $1\leq k\leq N$, the correlation matrix is simply $$ C_{mn} = \sum_{k=1}^N U^{*}_{km} U_{kn} $$.
Now back to your question. The first step is to compute the coefficients $V_{rs}$, where $$ a_r(t) = \sum_{s} V_{rs}a_s. $$ (Note also that operators with no time labelled means the operator at $t=0$). This can be computed via a matrix exponential $$ V = e^{-i H^{f}_{\mathrm{mat}} t}, $$ where $H^{f}_{\mathrm{mat}}$ is simply a tridiagonal matrix with off-diagonal elements $-J$. (This is exactly the first-quantized form of $H^{f}$). In your case, $V$ can be computed analytically by fully diagonalizing $H^{f}_{\mathrm{mat}}$ using Fourier modes. More generally, If $J$ also varies with position, then the matrix exponential must be computed numerically.
Then the correlation function is simply $$\langle a_r(t) a^{\dagger}_s(0)\rangle = \sum_{s'} V_{rs'} \langle a_{s'}a^{\dagger}_s\rangle = \sum_{s'} V_{rs'} (C_{ss'}+\delta_{ss'}),$$ where in the last equality I have used the canonical commutation relation to switch the creation and annihilation operators.