Multiplicity functions of quasars: Predictions from the MassiveBlackII simulation

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.