$\mathbf{Note :}$ Here equation of motion method is used to deduce the result OP is after.
With the OP's Hamiltonian:
$$\hat H = \sum_{i,j=1,1}^{N,N}\mathbb{H}_{ij}^{}\hat c_{i}^{\dagger} \hat c_{j}^{}$$
and the definition of OP's retarded Green's function as :
$$\mathbb{G}_{ij}^{}(t_{1}^{},t_{2}^{})=-\frac{i}{\hbar}\Theta(t_{1}^{}-t_{2}^{})\langle \{ \hat c_{i}^{}(t_{1}^{}),\hat c_{j}^{\dagger} (t_{2}^{})\} \rangle$$
assumed here is the standard variant, presuming OP's question has a typo.
Here $$\hat O(t)= e_{}^{\frac{i}{\hbar}\hat H t} \hat O e_{}^{-\frac{i}{\hbar}\hat H t}.$$
It can be shown that $\mathbb{G}_{}^{}(t_{1}^{},t_{2}^{})$ (matrix of $\mathbb{G}_{ij}^{}(t_{1}^{},t_{2}^{})$'s) satisfies the following differential equation's:
$$\left[i\hbar\frac{\partial}{\partial t_{1}^{}}\mathbb{I} -\mathbb{H}\right]\mathbb{G}_{}^{}(t_{1}^{},t_{2}^{})=\delta(t_{1}^{}-t_{2}^{})\mathbb{I}$$
and
$$\mathbb{G}_{}^{}(t_{1}^{},t_{2}^{})\left[-i\hbar\frac{\partial}{\partial t_{2}^{}}\mathbb{I} -\mathbb{H}\right]=\delta(t_{1}^{}-t_{2}^{})\mathbb{I}$$
here $\mathbb{I}$ is an identity matrix of dimension, $N \times N$.
The solution of above two equations should clearly be of the form (the infamous jump discontinuity) :
$$\mathbb{G}_{}^{}(t_{1}^{},t_{2}^{})=\Theta(t_{1}^{}-t_{2}^{})\mathbb{G}_{}^{>}(t_{1}^{},t_{2}^{})+\Theta(t_{2}^{}-t_{1}^{})\mathbb{G}_{}^{<}(t_{1}^{},t_{2}^{})$$
with (the magnitude of jump discontinuity)
$$\mathbb{G}_{}^{>}(t_{}^{},t_{}^{})-\mathbb{G}_{}^{<}(t_{}^{},t_{}^{}) = -\frac{i}{\hbar}\mathbb{I}$$
which can be deduced by integrating either first differential equation with respect to $t_1$ between $t_2+0^+$ and $t_2-0^+$ and using the above ansatz proposed followed by the identification $t_2=t$ or integrating second differential equation with respect to $t_2$ between $t_1+0^+$ and $t_1-0^+$ and using the above ansatz proposed followed by the identification $t_1=t$.
Further from the definition of $\mathbb{G}_{ij}^{}(t_{1}^{},t_{2}^{})$ it can be seen that $\mathbb{G}_{}^{<}(t_{1}^{},t_{2}^{})=\mathbb{O}$.
Hence $\mathbb{G}_{}^{>}(t_{}^{},t_{}^{}) = -\frac{i}{\hbar}\mathbb{I}$.
Using thus deduced results and the resultant ansatz i.e.,
$$\mathbb{G}_{}^{}(t_{1}^{},t_{2}^{})=\Theta(t_{1}^{}-t_{2}^{})\mathbb{G}_{}^{>}(t_{1}^{},t_{2}^{})$$
the following equations results,
$$\left[i\hbar\frac{\partial}{\partial t_{1}^{}}\mathbb{I} -\mathbb{H}\right]\mathbb{G}_{}^{>}(t_{1}^{},t_{2}^{})=\mathbb{O}$$
and
$$\mathbb{G}_{}^{>}(t_{1}^{},t_{2}^{})\left[-i\hbar\frac{\partial}{\partial t_{2}^{}}\mathbb{I} -\mathbb{H}\right]=\mathbb{O}$$
which can be easily integrated from $t$ to $t_{1}^{}$ and $t$ to $t_{2}^{}$ respectively to get
$$\mathbb{G}_{}^{>}(t_{1}^{},t_{2}^{})=e_{}^{-\frac{i}{h}\mathbb{H} (t_{1}^{}-t)}\mathbb{G}_{}^{>}(t_{}^{},t_{}^{})e_{}^{\frac{i}{h}\mathbb{H} (t_{2}^{}-t)}$$
which upon using $\mathbb{G}_{}^{>}(t_{}^{},t_{}^{}) = -\frac{i}{\hbar}\mathbb{I}$ gives:
$$\mathbb{G}_{}^{>}(t_{1}^{},t_{2}^{})=-\frac{i}{\hbar}e_{}^{-\frac{i}{h}\mathbb{H} (t_{1}^{}-t_{2}^{})}$$
using which finally the following falls out :
$$\mathbb{G}_{}^{}(t_{1}^{},t_{2}^{})=-\frac{i}{\hbar}\Theta(t_{1}^{}-t_{2}^{})e_{}^{-\frac{i}{h}\mathbb{H} (t_{1}^{}-t_{2}^{})}.$$
clearly $\mathbb{G}_{}^{}(t_{1}^{},t_{2}^{})$ is time translationally invariant, hence without loss off generality, the following can be considered instead:
$$\mathbb{G}_{}^{}(t)=-\frac{i}{\hbar}\Theta(t)e_{}^{-\frac{i}{h}\mathbb{H} t}.$$
Fourier transforming this as $\tilde{\mathbb{G}}_{}^{}(E)=\int_{-\infty}^{+\infty}dt\mathbb{G}_{}^{}(t)e_{}^{\frac{i}{\hbar}E t}$ finally reveals the sought after result :
$$\tilde{\mathbb{G}}_{}^{}(E)=\left[E\mathbb{I}-\mathbb{H}+i 0_{}^{+}\mathbb{I}\right]_{}^{-1}.$$
