Binary neutron star populations across cosmic time: the impact of binary stellar evolution uncertainties

We investigate the impact of uncertainties in binary stellar evolution on the formation efficiency and cosmological merger rate of BNS systems. In particular, we aim to determine whether robust signatures of binary evolution can be identified across a broad range of metallicities. We perform a systematic exploration of 72 binary evolution models using the population synthesis code COSMIC. We vary four key parameters controlling the onset and efficiency of common-envelope evolution, the natal kick prescription, and the treatment of kicks during the common envelope phase. For each model, we compute the BNS formation efficiency and cosmological merger rate using several prescriptions for the cosmic metallicity-redshift relation. We identify characteristic features in the metallicity dependence of the BNS formation efficiency that arise from specific evolutionary processes. These include the onset of first-giant-branch skipping at low metallicity and a non-monotonic evolution between Z=10-3 and 7.5 times 10-3. All models also predict a bimodal distribution of BNS progenitor masses associated with two distinct formation channels. The common envelope ejection efficiency and natal kick prescription remain the dominant sources of uncertainty, producing variations of several orders of magnitude in both the BNS formation efficiency and the cosmological merger rate. Finally, we show that different combinations of binary-evolution parameters can produce similar merger-rate histories, highlighting significant degeneracies among population-synthesis models. Current gravitational-wave observations can rule out only the least efficient formation scenarios. Breaking these degeneracies will require larger samples of BNS mergers expected from future third-generation gravitational-wave observatories, and improved observational constraints on the evolution of metallicity across cosmic time.