  All modern galaxy formation models employ stochastic elements in their
sub-grid prescriptions to discretise continuous equations across the time
domain. In this paper, we investigate how the stochastic nature of these
models, notably star formation, black hole accretion, and their associated
feedback, that act on small ($<$ kpc) scales, can back-react on macroscopic
galaxy properties (e.g. stellar mass and size) across long ($>$ Gyr)
timescales. We find that the scatter in scaling relations predicted by the
EAGLE model implemented in the SWIFT code can be significantly impacted by
random variability between re-simulations of the same object, even when
galaxies are resolved by tens of thousands of particles. We then illustrate how
re-simulations of the same object can be used to better understand the
underlying model, by showing how correlations between galaxy stellar mass and
black hole mass disappear at the highest black hole masses ($M_{\rm BH} > 10^8$
M$_\odot$), indicating that the feedback cycle may be interrupted by external
processes. We find that although properties that are collected cumulatively
over many objects are relatively robust against random variability (e.g. the
median of a scaling relation), the properties of individual galaxies (such as
galaxy stellar mass) can vary by up to 25\%, even far into the well-resolved
regime, driven by bursty physics (black hole feedback) and mergers between
galaxies. We suggest that studies of individual objects within cosmological
simulations be treated with caution, and that any studies aiming to closely
investigate such objects must account for random variability within their
results.
