  We examine multiple AGN systems (triples and quadruples, in particular) in
the \texttt{MassiveBlackII} simulation over a redshift range of $0.06\lesssim z
\lesssim 4$. We identify AGN systems (with bolometric luminosity
$L_{\mathrm{bol}}>10^{42}~\mathrm{ergs/sec}$) at different scales~(defined by
the maximum distance between member AGNs) to determine the AGN multiplicity
functions. This is defined as the volume/ surface density of AGN systems per
unit \textit{richness} $R$, the number of AGNs in a system. We find that
gravitationally bound multiple AGN systems tend to populate scales of
$\lesssim0.7~\mathrm{cMpc}/h$; this corresponds to angular separations of
$\lesssim100~\mathrm{arcsec}$ and a line of sight velocity difference
$\lesssim200~\mathrm{km/sec}$. The simulation contains $\sim 10$ and $\sim100$
triples/quadruples per $\mathrm{deg}^2$ up to depths of DESI ($g\lesssim24$)
and LSST ($g\lesssim26$) imaging respectively; at least $20\%$ of these should
be detectable in spectroscopic surveys. The simulated quasar
($L_{\mathrm{bol}}>10^{44}~\mathrm{ergs/sec}$) triples and quadruples
predominantly exist at $1.5\lesssim z \lesssim 3$. Their members have black
hole masses $10^{6.5}\lesssim M_{bh}\lesssim 10^{9}~M_{\odot}/h$ and live in
separate (one central and multiple satellite) galaxies with stellar masses
$10^{10}\lesssim M_{*}\lesssim 10^{12}~M_{\odot}/h$. They live in the most
massive haloes (for e.g. $\sim 10^{13}~M_{\odot}/h$ at $z=2.5$; $\sim
10^{14}~M_{\odot}/h$ at $z=1$) in the simulation. Their detections provide an
exciting prospect for understanding massive black hole growth and their merger
rates in galaxies in the era of multi-messenger astronomy.

Astrophysics of Galaxies;Cosmology and Nongalactic Astrophysics