The Monte Carlo method is widely used to compute the fundamental eigenfunction and eigenvalue for nuclear systems. However, the standard power iteration method computes only the fundamental eigenmode, while it would be beneficial to also compute the higher eigenfunctions and eigenvalues to support the reactor transient analysis, stability analysis, and assessments of nuclear safety, as well as to enable certain source convergence acceleration techniques. Modifications to the power method have been developed that in principle can accomplish this goal, but they typically lead to unphysical positive and negative particles requiring a procedure to compute the net-weight deposition. In this paper, we present a new mechanism that enables the Monte Carlo procedure, with certain modifications, to compute the second eigenfunction and eigenvalue for one-dimensional (1-D) problems. The method could in principle be extended to higher harmonics and general geometries. The results from numerical examples, including a 1-D, two-group, multiregion example, are consistent with reference results. Moreover, the extra computational cost of this method is of the same order of magnitude as the conventional Monte Carlo simulations. This method can be applied solely to solve for the high eigenmodes, or implemented as a part of a net-weight computation mechanism when negative particles are present in the modified power iteration method.