Inversion of Self-Potential (SP) data poses a significant challenge in exploration geophysics due to the inherent non-linearity and presence of multiple local optima. This study presents the first application of Harris Hawks Optimization (HHO) algorithm for quantitative interpretation of SP data based on a dipole model characterized by polarization intensity ( \(k\) ) and pole positions ( \(x_{1}\) , \(z_{1}\) , \(x_{2}\) , \(z_{2}\) ). Unlike conventional optimization methods, HHO’s unique two-phase mechanism—dynamically transitioning between exploration and exploitation through soft and hard besiege strategies—provides superior capability in escaping local optima, a critical advantage for the highly non-linear SP inverse problem. Validation on synthetic data (noise-free and 5–25% Gaussian noise) demonstrated exceptional accuracy, with relative errors below 2.54% for noise-free cases and \({R}^{2}=0.9725\) at 15% noise level. Comprehensive benchmarking against three established algorithms (PSO, GA, WOA) over 30 independent runs revealed HHO’s competitive performance (mean RMS = 25.012 mV, SD = 0.029 mV), achieving stability comparable to PSO while requiring no algorithm-specific parameter tuning. Parameter sensitivity analysis identified \({x}_{1}\) as the most influential parameter (normalized sensitivity = 1.0). Field validation on SP data from Surda copper mine, India, yielded excellent agreement (RMS = 5.69 mV, \({R}^{2}=0.9982\) ), with interpretations consistent with regional geology and previous studies. These results establish HHO as a robust, accurate, and practically viable tool for SP data inversion in mineral exploration.