  The forthcoming Laser Interferometer Space Antenna (LISA) will probe the
population of coalescing massive black hole (MBH) binaries up to the onset of
structure formation. Here we simulate the galactic-scale pairing of $\sim10^6
M_\odot$ MBHs in a typical, non-clumpy main-sequence galaxy embedded in a
cosmological environment at $z = 7-6$. In order to increase our statistical
sample, we adopt a strategy that allows us to follow the evolution of six
secondary MBHs concomitantly. We find that the magnitude of the
dynamical-friction induced torques is significantly smaller than that of the
large-scale, stochastic gravitational torques arising from the perturbed and
morphologically evolving galactic disc, suggesting that the standard dynamical
friction treatment is inadequate for realistic galaxies at high redshift. The
dynamical evolution of MBHs is very stochastic, and a variation in the initial
orbital phase can lead to a drastically different time-scale for the inspiral.
Most remarkably, the development of a galactic bar in the host system either
significantly accelerates the inspiral by dragging a secondary MBH into the
centre, or ultimately hinders the orbital decay by scattering the MBH in the
galaxy outskirts. The latter occurs more rarely, suggesting that galactic bars
overall promote MBH inspiral and binary coalescence. The orbital decay time can
be an order of magnitude shorter than what would be predicted relying on
dynamical friction alone. The stochasticity, and the important role of global
torques, have crucial implications for the rates of MBH coalescences in the
early Universe: both have to be accounted for when making predictions for the
upcoming LISA observatory.
