The scatter (σ sSFR) of the specific star formation rates of galaxies is a measure of the diversity in their star formation histories (SFHs) at a given mass. In this paper, we employ the Evolution and Assembly of GaLaxies and their Environments (EAGLE) simulations to study the dependence of the σ sSFR of galaxies on stellar mass (M ⋆) through the σ sSFR-M ⋆ relation in z ∼ 0-4. We find that the relation evolves with time, with the dispersion depending on both stellar mass and redshift. The models point to an evolving U-shaped form for the σ sSFR-M ⋆ relation, with the scatter being minimal at a characteristic mass M ⋆ of 109.5 M ☉ and increasing both at lower and higher masses. This implies that the diversity of SFHs increases toward both the low- and high-mass ends. We find that feedback from active galactic nuclei is important for increasing the σ sSFR for high-mass objects. On the other hand, we suggest that feedback from supernovae increases the σ sSFR of galaxies at the low-mass end. We also find that excluding galaxies that have experienced recent mergers does not significantly affect the σ sSFR-M ⋆ relation. Furthermore, we employ the EAGLE simulations in combination with the radiative transfer code SKIRT to evaluate the effect of SFR/stellar mass diagnostics in the σ sSFR-M ⋆ relation, and find that the SFR/M ⋆ methodologies (e.g., SED fitting, UV+IR, UV+IRX-β) widely used in the literature to obtain intrinsic properties of galaxies have a large effect on the derived shape and normalization of the σ sSFR-M ⋆ relation.