The evolution of galaxies is governed by equations with chaotic solutions: gravity and compressible (magneto-)hydrodynamics. While this microscale chaos and stochasticity has been well studied, it is poorly understood how it couples to macroscale properties examined in simulations of galaxy formation. In this paper, we use tiny perturbations introduced by floating-point round-off, random number generators, and seemingly trivial differences in algorithmic behaviour as seeds for chaotic behaviour. These can ultimately grow to produce non-trivial differences in star formation histories, circumgalactic medium properties, and the distribution of stellar mass. We examine the importance of stochasticity due to discreteness noise, variations in merger timings, and how self-regulation moderates the effects of this stochasticity. We show that chaotic variations in stellar mass can be suppressed through gas exhaustion, or be maintained at a roughly constant variance through stellar feedback. We also find that galaxy mergers are critical points from which large (as much as a factor of 2) variations can grow in global quantities such as the total stellar mass. These variations can grow and persist for Gyr before regressing towards the mean. These results show that detailed comparisons of simulations require serious consideration of the magnitude of effects compared to run-to-run chaotic variation, and may significantly complicate interpreting the impact of different physical models. Understanding the results of simulations require us to understand that the process of simulation is not a mapping of an infinitesimal point in configuration space to another, final infinitesimal point. Instead, simulations map a point in a space of possible initial conditions points to a volume of possible final states.