Stochastic partial differential equations (SPDEs) are often difficult to solve numerically due to their low regularity and high dimensionality. These challenges limit the practical use of computer-aided studies and pose significant barriers to statistical analysis of SPDEs. In this work, we introduce a highly efficient multi-index Monte Carlo method (MIMC) designed to approximate statistics of mild solutions to semilinear parabolic SPDEs. Key to our approach is the proof of a multiplicative convergence property for coupled solutions generated by an exponential integrator numerical solver, which we incorporate with MIMC. We further describe theoretically how the asymptotic computational cost of MIMC can be bounded in terms of the input accuracy tolerance, as the tolerance goes to zero. Notably, our methodology illustrates that for an SPDE with low regularity, MIMC offers substantial performance improvements over other viable methods. Numerical experiments comparing the performance of MIMC with the multilevel Monte Carlo method on relevant test problems validate our theoretical findings. These results also demonstrate that MIMC significantly outperforms state-of-the-art multilevel Monte Carlo, thereby underscoring its potential as a robust and tractable tool for solving semilinear parabolic SPDEs.