In statistical seismology, the Epidemic Type Aftershocks Sequence (ETAS) model is a branching process used world-wide to forecast earthquake intensity rates and reproduce many statistical features observed in seismicity catalogs. In this paper, we describe a fractional differential equation that governs the earthquake intensity rate of the pure temporal ETAS model by using the Caputo fractional derivative and we solve it analytically. We highlight that the tools and special functions of fractional calculus simplify the classical methods employed to obtain the intensity rate and let us describe the change of solution decay for large times. We also apply and discuss the theoretical results to the Japanese catalog in the period 1965-2003.