Abstract
We combine the swap Monte Carlo algorithm to long multi-CPU molecular dynamics simulations to analyse the equilibrium relaxation dynamics of model supercooled liquids over a time window covering ten orders of magnitude for temperatures down to the experimental glass transition temperature $T_g$. The analysis of \rev{several} time correlation functions coupled to spatio-temporal resolution of particle motion allow us to elucidate the nature of the equilibrium dynamics in deeply supercooled liquids. We find that structural relaxation starts at early times in rare localised regions characterised by a waiting time distribution that develops a power law near $T_g$. At longer times, relaxation events accumulate with increasing probability in these regions as $T_g$ is approached. This accumulation leads to a power-law growth of the linear extension of relaxed domains with time with a large, temperature-dependent dynamic exponent. Past the average relaxation time, unrelaxed domains slowly shrink with time due to relaxation events happening at their boundaries. Our results provide a complete microscopic description of the particle motion responsible for key experimental signatures of glassy dynamics, from the shape and temperature evolution of relaxation spectra to the core features of dynamic heterogeneity. They also provide a microscopic basis to understand the emergence of dynamic facilitation in deeply supercooled liquids and allow us to critically reassess theoretical descriptions of the glass transition.