What are the ways to carry out time propagation numerically? There are many ways to solve the time-dependent Schrodinger equation (TDSE) and find the wavefunction $|\Psi(t)\rangle$ for a given Hamiltonian. For example, consider a tight binding type Hamiltonian:
$$
\hat{H}
= \varepsilon \sum_{i}^{}{|i\rangle\langle i|}
+ \sum_{\left<i,j\right>}^{}{V_{ij}|i\rangle\langle j|}
$$
where $\{|i\rangle \}$ is the site basis.
To find the transition amplitude $\langle p | e^{-i\frac{H t}{\hbar}} | q \rangle$ (where $\{p,q\}$ belong to the site basis) which method is more efficient numerically and why?
One way is to find the eigenvalues and eigenvectors through exact diagonalization and further expand the wave function in terms of the eigenkets. The other way could be to do matrix exponentiation for the Hamiltonian matrix. What would be the case in case of other propagation methods like split operator method and Chebyshev operator method?
 A: I am not sure exactly what you are asking, but I can give a bit of detail about how we might implement a split operator method. I apologize if this is a orthogonal to exactly what you were asking. 
As you well know, time evolution in QM obeys the Schrodinger equation 
\begin{equation}
    i \hbar \frac{\partial\psi(t)}{\partial t} = \mathcal{H} \psi(t).
\end{equation}
Here, $\mathcal{H}$ is the Hamiltonian of the system in question. Usually, we can write $\mathcal{H}$ as the sum of a kinetic term and a few others, such as an interaction term or a external potential. Here I will focus only on the case of an external potential. In this case we may write
\begin{equation}
    \mathcal{H} = \frac{\bar{P}^2}{2m} + V(x, t)
\end{equation}
Time evolution in QM is unitary, so there is some unitary transformation $\hat{U}$ that takes $\psi(x, 0) \rightarrow \psi(x,t)$. We can plug this into the Schrodinger equation to see that:
\begin{align*}
        i \hbar \frac{\partial\psi(t)}{\partial t} &= \mathcal{H} \psi(t) \\
        i \hbar  \frac{\partial\hat{U} \psi(t)}{\partial t} &= \mathcal{H} \hat{U} \psi{(x, 0)}
\end{align*}
Since this is true for an arbitrary $\psi$, it must be true for $\hat{U}$ itself. Setting $\hbar=1$ for ease of notation, we then have the following relation:
\begin{equation}\label{eqn:unitschro}
    i \frac{\partial \hat{U}}{\partial t} = \mathcal{H} \hat{U}
\end{equation}
with solution
\begin{equation}
    \hat{U} = e^{-i \int_{0}^t \mathcal{H}(t') dt'}
\end{equation}
Thus if we initialize a state to $\psi(x,0)$ we may evolve it by multiplying by the matrix exponential above, or:
\begin{equation}
    \psi(x, t) = e^{-i \mathcal{H} t} \psi(x, 0).
\end{equation}
if $\mathcal{H}$ does not depend on $t$.
Now comes actually splitting the operator. Recall that we can write \begin{equation}
    \mathcal{H} = \frac{\bar{P}^2}{2m} + V(x, t)
\end{equation}
Recalling that $P \rightarrow -i\nabla$ in the position basis, and that in general $P$ and $V$ will not commute, we may write
\begin{equation}
    e^{-i \left(\frac{\bar{P}^2}{2m} + V(x, t)\right) \Delta t} = e^{\frac{i \Delta t}{4m}\nabla^2} e^{-i V \Delta t} e^{\frac{i \Delta t}{4m}\nabla^2} + \mathcal{O}(t^3) \label{eqn:splitoperator} 
\end{equation}
where $\Delta t$ is our discretised time step.
The higher order terms on the right come from the fact that the kinetic and potential terms will not commute. This can be seen by expanding the left hand side using the series definition of the exponential, and rewriting the result as a Dyson series. I believe it can also be shown using the Baker-Campell-Hausdorff theorem. In any case, the point is that our original matrix exponential can be written in the above form. This is useful for numerical evolving the wavefunction due to the following observation: In momentum space, the $e^{\frac{i\Delta t}{4m}\nabla^2}$ term is a diagonal matrix which is computationally inexpensive to multiply by, and in position space the same is true for $e^{-iV\Delta t}$. So a rough algorithm to evolve the initial state $\psi(x, 0)$ would be:
Initialize grid    
Initialize(psi(x,0))
Initialize (V(x,0))
psi(k) = FFT(\psi(x))
for each timestep:
    psi(k) = e^{\frac{i\ Delta t}{4m}k^2} * psi(k)
    psi(x) = InverseFFT(psi(k))
    psi(x) = e^{-iV\Delta t} * psi(x)
    psi(k) = FFT(psi(x))
    psi(k) =  e^{\frac{i\Delta t}{4m}k^2} * psi(k)
psi(x) = InverseFFT(psi(k))

Essentially, each time step consists of propagating the free particle of the same mass over half a time interval, multiplying by a phase shift that comes from the potential applied over the entire interval, and finishing the step off with a final half time step of free particle propagation. When you have multiple time steps strung together, the final half step of the previous step combines with the first half step of the current one, so you essentially just apply successive full-step propagations alternating between real and momentum space. In general this is not so bad to do since it is limited (usually) by how fact you can do an FFT for whatever grid size you are working with. Very efficient algorithms (FFTW for example) exist for doing this.
