  We model here the merger histories of the supermassive black hole (SMBH)
population in the late stages of a cosmological simulation of a $\sim 2 \times
10^{13} M_\odot$ galaxy group. The gravitational dynamics around the several
tens of SMBHs ($M_{\bullet} \gtrsim 7.5\times 10^7 M_\odot$) hosted by the
galaxies in the group is computed at high accuracy using regularized
integration with the KETJU code. The 11 SMBHs which form binaries and
hierarchical triplets eventually merge after hardening through dynamical
friction, stellar scattering, and gravitational wave (GW) emission. The
binaries form at eccentricities of $e \sim 0.3$-$0.9$, with one system evolving
to a very high eccentricity of $e = 0.998$, and merge on timescales of a few
tens to several hundred megayears. During the simulation the merger-induced GW
recoil kicks eject one SMBH remnant from the central host galaxy. This
temporarily drives the galaxy off the $M_{\bullet}$-$\sigma_{\star}$ relation,
however the galaxy returns to the relation due to subsequent galaxy mergers,
which bring in new SMBHs. This showcases a possible mechanism contributing to
the observed scatter of the $M_{\bullet}$-$\sigma_{\star}$ relation. Finally,
we show that Pulsar Timing Arrays and LISA would be able to detect parts of the
GW signals from the SMBH mergers that occur during the $\sim 4\,\mathrm{Gyr}$
time span simulated with KETJU.
