  A recently developed self-healing diffusion Monte Carlo algorithm
  [PRB 79, 195117] is extended to the calculation of excited states. The
formalism is based on an excited-state fixed-node approximation and the mixed
estimator of the excited-state probability density. The fixed-node ground state
wave-functions of inequivalent nodal pockets are found simultaneously using a
recursive approach. The decay of the wave-function into lower energy states is
prevented using two methods: i) The projection of the improved trial-wave
function into previously calculated eigenstates is removed. ii) The reference
energy for each nodal pocket is adjusted in order to create a kink in the
global fixed-node wave-function which, when locally smoothed out, increases the
volume of the higher energy pockets at the expense of the lower energy ones
until the energies of every pocket become equal. This reference energy method
is designed to find nodal structures that are local minima for arbitrary
fluctuations of the nodes within a given nodal topology. We demonstrate in a
model system that the algorithm converges to many-body eigenstates in
bosonic-like and fermionic cases.
