  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.
