A time averaged evolution operator for excited state calculation

Question

I created this program for calculating excited states for a quantum system with 2 degrees of freedom without having to find the lower-energy eigenfunctions of \(\hat{H}\) first. The same method should also work for systems with more position coordinates. There is a quite long documentation pdf file included there.

https://github.com/TeemuIsojarvi/Excited_state_solver

This should be faster than the usual method of finding the ground state of square of shifted Hamiltonian, \((\hat{H} - E_t)^2\), converging to the energy eigenvalue nearest to trial value \(E_t\). The squared Hamiltonian has fourth order spatial derivatives in position representation, which means that imaginary time evolution can't be done using only tridiagonal matrices even with operator splitting (unlike the numerical method in my application).

Here is the probability density function calculated for an antibonding orbital of a 2D hydrogen molecular ion

Has any numerical method like this been presented earlier? The time averaged evolution operator is mentioned here, but in different context: https://arxiv.org/pdf/2601.02977