Back-projection (BP) is a cornerstone method for imaging earthquake ruptures, particularly effective at teleseismic distances for deciphering large earthquake kinematics. Its superior resolution is attributed to the ability to resolve high-frequency (>1 Hz) seismic signals, where waveforms immediately following the first coherent arrivals are composed of waves scattered by small-scale seismic velocity heterogeneities. This scattering leads to waveform incoherence between neighboring stations, a phenomenon not captured by synthetic tests of BP using Green’s functions (GF) derived from oversimplified 1D or smooth 3D velocity models. Addressing this gap, we introduce a novel approach to generate synthetic Incoherent Green’s Functions (IGF) that include scattered waves, accurately mimicking the observed inter-station waveform coherence decay spatially and temporally. Our methodology employs a waveform simulator that adheres to ray theory for the travel times of scattered waves, aggregating them as incident plane waves to simulate the high-frequency scattered wavefield across a seismic array. Contrary to conventional views that scattered waves degrade BP imaging quality by reducing array coherence, our synthetic tests reveal that IGFs are indispensable for accurately imaging extensive ruptures. Specifically, the rapid decay of IGF coherence prevents early rupture segments from overshadowing subsequent ones, a critical flaw when using coherent GFs. By leveraging IGFs, we delve into previously unexplored aspects of BP imaging’s resolvability, sensitivity, fidelity, and uncertainty. Our investigation not only highlights and explains the commonly observed “tailing” and “shadowing” artefacts but also proposes a robust framework for identifying different rupture stages and quantifying their uncertainties, thereby significantly enhancing BP imaging accuracy.