The association of GRB170817A with a binary neutron star (BNS) merger has revealed that BNSs produce at least a fraction of short gamma-ray bursts (SGRBs). As gravitational wave (GW) detectors push their horizons, it is important to assess coupled electromagnetic (EM)/GW probabilities and maximize observational prospects. Here, we perform BNS population synthesis calculations with the code mobse, seeding the binaries in galaxies at three representative redshifts, z = 0.01, 0.1, and 1 of the Illustris TNG50 simulation. The binaries are evolved and their locations numerically tracked in the host galactic potentials until merger. Adopting the microphysics parameters of GRB170817A, we numerically compute the broad-band light curves of jets from BNS mergers, with the afterglow brightness dependent on the local medium density at the merger site. We perform Monte Carlo simulations of the resulting EM population assuming either a random viewing angle with respect to the jet, or a jet aligned with the orbital angular momentum of the binary, which biases the viewing angle probability for GW-triggered events. We find a gamma-ray detection probability of ∼ 2 per cent,10 per cent, and 40 per cent for BNSs at z = 1, 0.1, and 0.01, respectively, for the random case, rising to ∼ 75 per cent for the z = 0.01, GW-triggered aligned case. Afterglow detection probabilities of GW-triggered BNS mergers vary in the range of ∼ 0.3 - 0.5 per cent, with higher values for aligned jets, and are comparable across the high- and low-energy bands, unlike gamma-ray-triggered events (cosmological SGRBs) which are significantly brighter at higher energies. We further quantify observational biases with respect to host galaxy masses.
