  We report an empirical determination of the probability density functions
P(r) of the number r of earthquakes in finite space-time windows for the
California catalog, over fixed spatial boxes 5 x 5 km^2 and time intervals dt
=1, 10, 100 and 1000 days. We find a stable power law tail P(r) ~ 1/r^{1+mu}
with exponent mu \approx 1.6 for all time intervals. These observations are
explained by a simple stochastic branching process previously studied by many
authors, the ETAS (epidemic-type aftershock sequence) model which assumes that
each earthquake can trigger other earthquakes (``aftershocks''). An aftershock
sequence results in this model from the cascade of aftershocks of each past
earthquake. We develop the full theory in terms of generating functions for
describing the space-time organization of earthquake sequences and develop
several approximations to solve the equations. The calibration of the theory to
the empirical observations shows that it is essential to augment the ETAS model
by taking account of the pre-existing frozen heterogeneity of spontaneous
earthquake sources. This seems natural in view of the complex multi-scale
nature of fault networks, on which earthquakes nucleate. Our extended theory is
able to account for the empirical observation satisfactorily. In particular,
the adjustable parameters are determined by fitting the largest time window
$dt=1000$ days and are then used as frozen in the formulas for other time
scales, with very good agreement with the empirical data.
