We constructed a seismic source model for the 2015 MW 8.3 Illapel, Chile earthquake, which was carried out with the kinematic waveform inversion method adopting a novel inversion formulation that takes into account the uncertainty in the Green’s function, together with the hybrid backprojection method enabling us to track the spatiotemporal distribution of high-frequency (0.3–2.0 Hz) sources at high resolution by using globally observed teleseismic P-waveforms. A maximum slip amounted to 10.4 m in the shallow part of the seismic source region centered 72 km northwest of the epicenter and generated a following tsunami inundated along the coast. In a gross sense, the rupture front propagated almost unilaterally to northward from the hypocenter at \2 km/s, however, in detail the spatiotemporal slip distribution also showed a complex rupture propagation pattern: two up-dip rupture propagation episodes, and a secondary rupture episode may have been triggered by the strong high-frequency radiation event at the down-dip edge of the seismic source region. High-frequency sources tends to be distributed at deeper parts of the slip area, a pattern also documented in other subduction zone megathrust earthquakes that may reflect the heterogeneous distribution of fracture energy or stress drop along the fault. The weak excitation of high-frequency radiation at the termination of rupture may represent the gradual deceleration of rupture velocity at the transition zone of frictional property or stress state between the megathrust rupture zone and the swarm area.