  We discuss the latest results of numerical simulations following the orbital
decay of massive black hole pairs in galaxy mergers. We highlight important
differences between gas-poor and gas-rich hosts, and between orbital evolution
taking place at high redshift as opposed to low redshift. Two effects have a
huge impact and are rather novel in the context of massive black hole binaries.
The first is the increase in characteristic density of galactic nuclei of
merger remnants as galaxies are more compact at high redshift due to the way
dark halo collapse depends on redshift. This leads naturally to hardening
timescales due to 3-body encounters that should decrease by two orders of
magnitude up to $z=4$. It explains naturally the short binary coalescence
timescale, $\sim 10$ Myr, found in novel cosmological simulations that follow
binary evolution from galactic to milliparsec scales. The second one is the
inhomogeneity of the interstellar medium in massive gas-rich disks at high
redshift. In the latter star forming clumps 1-2 orders of magnitude more
massive than local Giant Molecular Clouds (GMCs) can scatter massive black
holes out of the disk plane via gravitational perturbations and direct
encounters. This renders the character of orbital decay inherently stochastic,
often increasing orbital decay timescales by as much as a Gyr. At low redshift
a similar regime is present at scales of $1-10$ pc inside Circumnuclear Gas
Disks (CNDs). In CNDs only massive black holes with masses below $10^7
M_{\odot}$ can be significantly perturbed. They decay to sub-pc separations in
up to $\sim 10^8$ yr rather than the in just a few million years as in a smooth
CND. Finally implications for building robust forecasts of LISA event rates are
discussed
