The detection of electromagnetic radiation (EM) accompanying the gravitational wave (GW) signal from the binary neutron star (BNS) merger GW170817 has revealed that these systems constitute at least a fraction of the progenitors of short gammaray bursts (SGRBs). As gravitational wave detectors keep pushing their detection horizons, it is important to assess coupled GW/EM probabilities, and how to maximize observational prospects. Here we perform population synthesis calculations of BNS evolution with the code MOBSE, and seed the binaries in galaxies at three representative redshifts (𝑧 = 0.01, 0.1, 1) of the Illustris TNG50 simulation. The binaries are evolved and their locations numerically tracked in the host galactic potentials until merger. Adopting the astrophysical parameters of GRB170817A as a prototype, we numerically compute the broadband lightcurves of jets from BNS mergers, with the afterglow brightness depending on the local medium density at the merger sites. 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 that ∼ 70 − 80% of BNSs from 𝑧 = 0.01 should be detectable in gamma-rays. The afterglow detection probabilities of GW-triggered BNS mergers vary between ∼ 0.3 − 0.7%, with higher values for jets aligned with the BNS angular momentum, and are comparable across the high and low-energy bands, unlike 𝛾-ray-triggered events (cosmological SGRBs) which are significantly brighter at higher energies. We further quantify observational biases with respect to host galaxy masses.