  The rate and location of Binary Neutron Star (BNS) mergers are determined by
a combination of the star formation history and the Delay Time Distribution
(DTD) function. In this paper, we couple the star formation rate histories
(SFRHs) from the IllustrisTNG model to a series of varied assumptions for the
BNS DTD to make predictions for the BNS merger host galaxy mass function. These
predictions offer two outcomes: (i) in the near term: influence BNS merger
event follow-up strategy by scrutinizing where most BNS merger events are
expected to occur and (ii) in the long term: constrain the DTD for BNS merger
events once the host galaxy mass function is observationally well determined.
From our fiducial model analysis, we predict that 50% of BNS mergers will occur
in host galaxies with stellar mass between $10^{10} - 10^{11}$ $M_\odot$, 68%
between $4 \times 10^{9} - 3\times 10^{11}$ $M_\odot$, and 95% between $4
\times 10^8 - 2 \times 10^{12}$ $M_\odot$. We find that the details of the DTD
employed does not have a strong effect on the peak of the host mass function.
However, varying the DTD provides enough spread that the true DTD can be
determined from enough electromagnetic observations of BNS mergers. Knowing the
true DTD can help us determine the prevalence of BNS systems formed through
highly eccentric and short separation fast-merging channels and can constrain
the dominant source of r-process material.
