  Massive black hole (MBH) binaries can form following a galaxy merger, but
this may not always lead to a MBH binary merger within a Hubble time. The
merger timescale depends on how efficiently the MBHs lose orbital energy to the
gas and stellar background, and to gravitational waves (GWs). In systems where
these mechanisms are inefficient, the binary inspiral time can be long enough
for a subsequent galaxy merger to bring a third MBH into the system. In this
work, we identify and characterize the population of triple MBH systems in the
Illustris cosmological hydrodynamic simulation. We find a substantial
occurrence rate of triple MBH systems: in our fiducial model, 22% of all binary
systems form triples, and $>70$% of these involve binaries that would not
otherwise merge by $z=0$. Furthermore, a significant subset of triples (6% of
all binaries, or more than a quarter of all triples) form a triple system at
parsec scales, where the three BHs are most likely to undergo a strong
three-body interaction. Crucially, we find that the rate of triple occurrence
has only a weak dependence on key parameters of the binary inspiral model
(binary eccentricity and stellar loss-cone refilling rate). We also do not
observe strong trends in the host galaxy properties for binary versus triple
MBH populations. Our results demonstrate the potential for triple systems to
increase MBH merger rates, thereby enhancing the low-frequency GW signals
detectable with pulsar timing arrays and with LISA.
