  Massive black hole binaries (MBHBs) are expected to form at the centre of
merging galaxies during the hierarchical assembly of the cosmic structure, and
are expected to be the loudest sources of gravitational waves (GWs) in the low
frequency domain. However, because of the dearth of energy exchanges with
background stars and gas, many of these MBHBs may stall at separations too
large for GW emission to drive them to coalescence in less than a Hubble time.
Triple MBH systems are then bound to form after a further galaxy merger,
triggering a complex and rich dynamics that can eventually lead to MBH
coalescence. Here we report on the results of a large set of numerical
simulations, where MBH triplets are set in spherical stellar potentials and MBH
dynamics is followed through 2.5 post-Newtonian order in the equations of
motion. From our full suite of simulated systems we find that a fraction
$\simeq 20-30$% of the MBH binaries that would otherwise stall are led to
coalesce within a Hubble time. The corresponding coalescence timescale peaks
around 300 Myr, while the eccentricity close to the plunge, albeit small, is
non-negligible ($\lesssim 0.1$). We construct and discuss marginalised
probability distributions of the main parameters involved and, in a companion
paper of the series, we will use the results presented here to forecast the
contribution of MBH triplets to the GW signal in the nHz regime probed by
Pulsar Timing Array experiments.
