There are various ways to numerically find the ground state energy and wavefunction of a many-body Hamiltonian. You can diagonalize the Hamiltonian and pick out the lowest eigenstate, or you use Lancoz.
My proposal is that can I use imaginary time propagation for many-body problems?
Make it simple, I multiply a trial wavefunction $$|\psi_0\rangle = \sum_{i\in\text{eigen}} w^0_i |i\rangle$$ by an operator $\exp(-\hat H \tau)$. Then we have $$|\psi_\tau\rangle = \sum_{i\in\text{eigen}} w^\tau_i |i\rangle$$ with $w_i^\tau = w^0_i\exp(-E_i\tau)$.
For a sufficiently large $\tau$, we have $\exp(-\hat H \tau)\approx w_g^\tau \hat P_g$. The trial function will be projected to the ground state!
Choosing a complete set of basis states, we can numerically calculate the operator by taking Taylor expansion of the exponential and iterate to the $n$-th order, finally we obtain a matrix. Now multiply the matrix on a trial function written in terms of the basis we chose, then normalize it, and we get the ground state wavefunction.
Will it be accurate, stable and fast?