Scaling of dynamic rupture processes from small to large earthquakes is critical to seismic hazard assessment. Large subduction earthquakes are typically remote, and we mostly rely on teleseismic body waves to extract information on their slip rate functions. We estimate the P wave source spectra of 942 thrust earthquakes of magnitude Mw 5.5 and above by carefully removing wave propagation effects (geometrical spreading, attenuation, and free surface effects). The conventional spectral model of a single‐corner frequency and high‐frequency falloff rate does not explain our data, and we instead introduce a double‐corner‐frequency model, modified from the Haskell propagating source model, with an intermediate falloff of f−1. The first corner frequency f1 relates closely to the source duration T1, its scaling follows
M0∝T13 for Mw<7.5, and changes to
M0∝T12 for larger earthquakes. An elliptical rupture geometry better explains the observed scaling than circular crack models. The second time scale T2 varies more weakly with moment,
M0∝T25, varies weakly with depth, and can be interpreted either as expressions of starting and stopping phases, as a pulse‐like rupture, or a dynamic weakening process. Estimated stress drops and scaled energy (ratio of radiated energy over seismic moment) are both invariant with seismic moment. However, the observed earthquakes are not self‐similar because their source geometry and spectral shapes vary with earthquake size. We find and map global variations of these source parameters.