From the Milky Way to the Early Universe: How Metallicity Shapes the Formation of Double Neutron Stars

Author(s)

Chattaraj, Abhishek, Andrews, Jeff J., Gallegos-Garcia, Monica, Briel, Max M., Zevin, Michael, Fragos, Tassos, Zapartas, Emmanouil, Aguilera-Dena, David R., Gossage, Seth, Ray, Anarya, Srivastava, Philipp M., Xing, Zepei

Abstract

Among the compact binaries detected through gravitational wave (GW) emission, double neutron stars (DNSs) hold the unique distinction of having well-characterized local analogs in the Milky Way. However, understanding their broader cosmological population remains challenging due to the complex metallicity-dependent evolution of their progenitors. Using the detailed binary population synthesis code POSYDON, we model DNS formation across a cosmological range of metallicities. We find that while DNSs are predominantly formed through common envelope evolution, metallicity-dependent radial expansion and stellar wind mass loss bifurcates the pathway into distinct subchannels, producing systematically distinct merging populations at different metallicities. Notably, our best representative Galactic models predict local merger rates in broad agreement with the latest GWTC-5 constraints, suggesting that the populations observed via radio pulsar surveys and through GW emission need not be in tension. Contrary to previous studies, we demonstrate that the intrinsic merger efficiency increases toward low metallicities, a trend robust across varying binary evolution assumptions, although the magnitude and the underlying stellar processes driving this enhancement vary across metallicity regimes. When convolved with the cosmic star formation history, our models reveal that the majority of merging DNSs originate from progenitors in the range $\sim 0.1\,Z_\odot - Z_\odot$, with significant supersolar contributions. Furthermore, we show that weaker stellar winds at low metallicities produce DNS mergers with systematically more symmetric mass ratios ($q \gtrsim0.9$). Additionally, we assess detectability within proposed configurations for third-generation GW observatories. Our study provides a framework for interpreting the cosmological DNS population in the upcoming multi-messenger era.

Figures

The radial expansion of a $\simeq 10.2\Msun$ H-rich star at five representative metallicities. At lower metallicities, the star remains relatively compact during the Hertzsprung Gap (HG) phase, delaying significant expansion until advanced evolutionary stages. Consequently, donors in this regime are highly evolved at the onset of RLO. On the other hand, higher metallicity models experience rapid HG expansion immediately after terminal age main sequence (TAMS), allowing donors to initiate RLO early in their post-main-sequence evolution as well as during the more evolved giant branches.
Caption The radial expansion of a $\simeq 10.2\Msun$ H-rich star at five representative metallicities. At lower metallicities, the star remains relatively compact during the Hertzsprung Gap (HG) phase, delaying significant expansion until advanced evolutionary stages. Consequently, donors in this regime are highly evolved at the onset of RLO. On the other hand, higher metallicity models experience rapid HG expansion immediately after terminal age main sequence (TAMS), allowing donors to initiate RLO early in their post-main-sequence evolution as well as during the more evolved giant branches.
Evolution of the wind mass-loss rate (top panel) and effective temperature (bottom panel) for a $\simeq 10.2\,\Msun$ H-rich model across five initial metallicities. The horizontal black line in the bottom panel (and its inset) denotes the bi-stability jump at $T_{\mathrm{eff}} \approx 25,000\,\mathrm{K}$ (see Section \ref{subsubsec:wind_mass_loss}). The inset panels isolate the late-stage, post-main-sequence evolution of the $Z=10^{-2}\Zsun$ and $0.1\Zsun$ models. As high-metallicity stars ($Z \gtrsim 0.1\,\Zsun$) expand and cool on the main sequence, they cross the bi-stability jump and experience a surge in the wind mass loss rate. Furthermore, as highlighted in the top inset, there is a ${\sim}2\,\mathrm{dex}$ enhancement in the supergiant wind mass-loss rate for $Z \gtrsim 0.1\Zsun$ compared to metal-poor models, which remain hotter and more compact. Ultimately, this extreme mass loss forces significant orbital widening in these high-metallicity binaries prior to RLO, leading to a distinct outcropping of models in our \texttt{CO-HMS} grid (Figure \ref{fig:CO-HMS_RLO_multi_panel}), responsible for the bifurcation within the CE formation channel.
Caption Evolution of the wind mass-loss rate (top panel) and effective temperature (bottom panel) for a $\simeq 10.2\,\Msun$ H-rich model across five initial metallicities. The horizontal black line in the bottom panel (and its inset) denotes the bi-stability jump at $T_{\mathrm{eff}} \approx 25,000\,\mathrm{K}$ (see Section \ref{subsubsec:wind_mass_loss}). The inset panels isolate the late-stage, post-main-sequence evolution of the $Z=10^{-2}\Zsun$ and $0.1\Zsun$ models. As high-metallicity stars ($Z \gtrsim 0.1\,\Zsun$) expand and cool on the main sequence, they cross the bi-stability jump and experience a surge in the wind mass loss rate. Furthermore, as highlighted in the top inset, there is a ${\sim}2\,\mathrm{dex}$ enhancement in the supergiant wind mass-loss rate for $Z \gtrsim 0.1\Zsun$ compared to metal-poor models, which remain hotter and more compact. Ultimately, this extreme mass loss forces significant orbital widening in these high-metallicity binaries prior to RLO, leading to a distinct outcropping of models in our \texttt{CO-HMS} grid (Figure \ref{fig:CO-HMS_RLO_multi_panel}), responsible for the bifurcation within the CE formation channel.
Slices of the high resolution supplementary set of \texttt{CO-HMS} models terminating in unstable RLO across eight metallicities. The color coding denotes the donor's carbon-oxygen core mass at the onset of RLO. At low metallicities ($Z \lesssim 10^{-2}\,Z_\odot$), suppressed radial expansion during the HG phase ensures that all donors are highly evolved with developed C/O cores (Section \ref{subsubsec:radius_evol}). At higher metallicities ($Z \gtrsim 0.1\,Z_\odot$), two distinct physical effects drive a bifurcation in the parameter space. First, enhanced radial expansion allows less-evolved donors to fill their Roche lobes before developing a C/O core (dark purple models, $M_{\mathrm{CO-core}} = 0$), leading to DNS formation through the \Hecore formation channel. Second, extreme supergiant wind mass loss strips the envelopes of the highly evolved donors, forcing significant orbital widening prior to RLO (Section \ref{subsubsec:wind_mass_loss}). This supergiant wind-driven evolution shifts the corresponding \COcore channel into a distinct upper-left outcropping characterized by longer orbital periods and lower donor masses (Section \ref{subsec:CO-HMS_RLO}). Ultimately, the orbital expansion is so severe for $Z \gtrsim \Zsun$ (note the different y-axis ranges for these panels) that the \COcore channel fails to produce any BNSs that merge (Section \ref{subsec:all_form_channels}).
Caption Slices of the high resolution supplementary set of \texttt{CO-HMS} models terminating in unstable RLO across eight metallicities. The color coding denotes the donor's carbon-oxygen core mass at the onset of RLO. At low metallicities ($Z \lesssim 10^{-2}\,Z_\odot$), suppressed radial expansion during the HG phase ensures that all donors are highly evolved with developed C/O cores (Section \ref{subsubsec:radius_evol}). At higher metallicities ($Z \gtrsim 0.1\,Z_\odot$), two distinct physical effects drive a bifurcation in the parameter space. First, enhanced radial expansion allows less-evolved donors to fill their Roche lobes before developing a C/O core (dark purple models, $M_{\mathrm{CO-core}} = 0$), leading to DNS formation through the \Hecore formation channel. Second, extreme supergiant wind mass loss strips the envelopes of the highly evolved donors, forcing significant orbital widening prior to RLO (Section \ref{subsubsec:wind_mass_loss}). This supergiant wind-driven evolution shifts the corresponding \COcore channel into a distinct upper-left outcropping characterized by longer orbital periods and lower donor masses (Section \ref{subsec:CO-HMS_RLO}). Ultimately, the orbital expansion is so severe for $Z \gtrsim \Zsun$ (note the different y-axis ranges for these panels) that the \COcore channel fails to produce any BNSs that merge (Section \ref{subsec:all_form_channels}).
Slices of the \texttt{CO-HMS\_RLO} grids at $Z = 10^{-4}\,Z_\odot$ (left), $0.1\,Z_\odot$ (middle), and $Z_\odot$ (right), for a $1.43\,M_\odot$ NS and a H-rich donor at the onset of RLO. The grids combine the standard \posydon DR2 grids with our supplementary models (Section~\ref{subsec:setup_modifications}). Diamond markers denote systems that undergo CE evolution and square markers represent stable RLO. Contours show the 75th-percentile distributions of the \Hecore (solid) and \COcore (dashed) subchannels in \texttt{MODEL03}.
Caption Slices of the \texttt{CO-HMS\_RLO} grids at $Z = 10^{-4}\,Z_\odot$ (left), $0.1\,Z_\odot$ (middle), and $Z_\odot$ (right), for a $1.43\,M_\odot$ NS and a H-rich donor at the onset of RLO. The grids combine the standard \posydon DR2 grids with our supplementary models (Section~\ref{subsec:setup_modifications}). Diamond markers denote systems that undergo CE evolution and square markers represent stable RLO. Contours show the 75th-percentile distributions of the \Hecore (solid) and \COcore (dashed) subchannels in \texttt{MODEL03}.
The intrinsic BNS merger efficiency per unit stellar mass formed as a function of progenitor metallicities for the full set of population synthesis models considered in this work (Table \ref{tab:model_populations}). Individual \texttt{POSYDON} model variations exploring uncertainties in binary physics are shown in thin grey lines, and our observationally calibrated subset is highlighted in color. For comparison, the dashed red line denotes the meta-analysis from \citet{vanSon+2025}, which consists of predictions from rapid population synthesis studies and only covers $Z \ge 10^{-2}\Zsun$. Our observationally calibrated models predict a decline in the efficiency toward higher metallicities, unlike the flatter trend of the meta-analysis over the metallicity range in common.
Caption The intrinsic BNS merger efficiency per unit stellar mass formed as a function of progenitor metallicities for the full set of population synthesis models considered in this work (Table \ref{tab:model_populations}). Individual \texttt{POSYDON} model variations exploring uncertainties in binary physics are shown in thin grey lines, and our observationally calibrated subset is highlighted in color. For comparison, the dashed red line denotes the meta-analysis from \citet{vanSon+2025}, which consists of predictions from rapid population synthesis studies and only covers $Z \ge 10^{-2}\Zsun$. Our observationally calibrated models predict a decline in the efficiency toward higher metallicities, unlike the flatter trend of the meta-analysis over the metallicity range in common.
The relative BNS merger efficiency across different progenitor metallicity regimes for all model populations. The panels provide the ratio of merger efficiencies between $10^{-4}\,Z_\odot$ and $0.1\,Z_\odot$ (\textit{left}), $0.1\,Z_\odot$ and $Z_\odot$ (\textit{middle}), and across the entire metallicity range considered in this work from $10^{-4}\,Z_\odot$ to $Z_\odot$ (\textit{right}). The horizontal dashed line indicates a value of unity, where the formation efficiency is equal between the two compared metallicities; models falling above the line favors BNSs merging in the relatively lower metallicity environment. \textit{Left panel:} The scatter of all models above and below unity indicates that efficiencies in the very metal-poor regime are sensitive to binary physics uncertainties, though all observationally calibrated models favor $10^{-4}\,Z_\odot$. \textit{Middle panel:} In stark contrast, every model falls strictly above unity, demonstrating a universal preference for $0.1\,Z_\odot$ over $Z_\odot$. This is driven by the suppression of the \COcore channel in forming merging BNS toward high metallicity due to extreme supergiant wind mass loss and orbital widening, increasing the binary's susceptibility to natal kick disruption. \textit{Right panel:} Across the full metallicity range, almost all models favor metal-poor environments. Our observationally calibrated subset predicts an efficiency enhancement of roughly an order of magnitude at $10^{-4}\,Z_\odot$ compared to $Z_\odot$.
Caption The relative BNS merger efficiency across different progenitor metallicity regimes for all model populations. The panels provide the ratio of merger efficiencies between $10^{-4}\,Z_\odot$ and $0.1\,Z_\odot$ (\textit{left}), $0.1\,Z_\odot$ and $Z_\odot$ (\textit{middle}), and across the entire metallicity range considered in this work from $10^{-4}\,Z_\odot$ to $Z_\odot$ (\textit{right}). The horizontal dashed line indicates a value of unity, where the formation efficiency is equal between the two compared metallicities; models falling above the line favors BNSs merging in the relatively lower metallicity environment. \textit{Left panel:} The scatter of all models above and below unity indicates that efficiencies in the very metal-poor regime are sensitive to binary physics uncertainties, though all observationally calibrated models favor $10^{-4}\,Z_\odot$. \textit{Middle panel:} In stark contrast, every model falls strictly above unity, demonstrating a universal preference for $0.1\,Z_\odot$ over $Z_\odot$. This is driven by the suppression of the \COcore channel in forming merging BNS toward high metallicity due to extreme supergiant wind mass loss and orbital widening, increasing the binary's susceptibility to natal kick disruption. \textit{Right panel:} Across the full metallicity range, almost all models favor metal-poor environments. Our observationally calibrated subset predicts an efficiency enhancement of roughly an order of magnitude at $10^{-4}\,Z_\odot$ compared to $Z_\odot$.
: SFR: IllustrisTNG
Caption : SFR: IllustrisTNG
: SFR: Madau+Fragos17
Caption : SFR: Madau+Fragos17
The BNS merger rate density evolution for all the population synthesis models in our work assuming two cosmic star formation histories: \texttt{IllustrisTNG} (left) and \texttt{MF17} (right). The grey lines representing our full set of population models with varying assumptions of binary physics exhibit a broad range in the rate evolution. Nevertheless, all models closely track the shape of the underlying SFH (blue dotted line), with merger rates typically peaking at $z\simeq2-3$ corresponding to cosmic noon. Additionally, our observationally calibrated subset (colored lines) broadly matches the latest local rate constraints from GWTC-5 (shaded in red, with an extended range up to $z=0.5$ for visualization purposes), although they lie near the upper end of the inferred range of $5.1-154.7\,\rm{Gpc^{-3}}\,\rm{yr^{-1}}$. This overlap may indicate that the observed Galactic DNS population need not be in tension with the merging BNS population inferred from GW observations.
Caption The BNS merger rate density evolution for all the population synthesis models in our work assuming two cosmic star formation histories: \texttt{IllustrisTNG} (left) and \texttt{MF17} (right). The grey lines representing our full set of population models with varying assumptions of binary physics exhibit a broad range in the rate evolution. Nevertheless, all models closely track the shape of the underlying SFH (blue dotted line), with merger rates typically peaking at $z\simeq2-3$ corresponding to cosmic noon. Additionally, our observationally calibrated subset (colored lines) broadly matches the latest local rate constraints from GWTC-5 (shaded in red, with an extended range up to $z=0.5$ for visualization purposes), although they lie near the upper end of the inferred range of $5.1-154.7\,\rm{Gpc^{-3}}\,\rm{yr^{-1}}$. This overlap may indicate that the observed Galactic DNS population need not be in tension with the merging BNS population inferred from GW observations.
The BNS merger rate density evolution for all the population synthesis models in our work assuming two cosmic star formation histories: \texttt{IllustrisTNG} (left) and \texttt{MF17} (right). The grey lines representing our full set of population models with varying assumptions of binary physics exhibit a broad range in the rate evolution. Nevertheless, all models closely track the shape of the underlying SFH (blue dotted line), with merger rates typically peaking at $z\simeq2-3$ corresponding to cosmic noon. Additionally, our observationally calibrated subset (colored lines) broadly matches the latest local rate constraints from GWTC-5 (shaded in red, with an extended range up to $z=0.5$ for visualization purposes), although they lie near the upper end of the inferred range of $5.1-154.7\,\rm{Gpc^{-3}}\,\rm{yr^{-1}}$. This overlap may indicate that the observed Galactic DNS population need not be in tension with the merging BNS population inferred from GW observations.
Caption The BNS merger rate density evolution for all the population synthesis models in our work assuming two cosmic star formation histories: \texttt{IllustrisTNG} (left) and \texttt{MF17} (right). The grey lines representing our full set of population models with varying assumptions of binary physics exhibit a broad range in the rate evolution. Nevertheless, all models closely track the shape of the underlying SFH (blue dotted line), with merger rates typically peaking at $z\simeq2-3$ corresponding to cosmic noon. Additionally, our observationally calibrated subset (colored lines) broadly matches the latest local rate constraints from GWTC-5 (shaded in red, with an extended range up to $z=0.5$ for visualization purposes), although they lie near the upper end of the inferred range of $5.1-154.7\,\rm{Gpc^{-3}}\,\rm{yr^{-1}}$. This overlap may indicate that the observed Galactic DNS population need not be in tension with the merging BNS population inferred from GW observations.
The redshift evolution of the BNS merger rate density for one of our observationally calibrated population $\texttt{MODEL03}$, assuming the \texttt{IllustrisTNG} star formation history. The total rate (black) is broken down into contributions from eight birth metallicities and compared against the GWTC-5 constraints on the local rate (gray shaded region, which has been extended to $z=0.5$ for visualization). Three features of the cosmological BNS population can be inferred: (i) the local Universe ($z < 0.1$) is overwhelmingly dominated by progenitors born in the $0.1 - 2\,Z_\odot$ range; (ii) the contribution of metal-poor progenitors ($Z < 0.1\,Z_\odot$) peak at higher redshifts ($z\gtrsim3$); and finally, (iii) due to relatively short delay times, the the merger rate evolution broadly follows the cosmic star formation history, making BNS mergers potential tracers of cosmic chemical evolution.
Caption The redshift evolution of the BNS merger rate density for one of our observationally calibrated population $\texttt{MODEL03}$, assuming the \texttt{IllustrisTNG} star formation history. The total rate (black) is broken down into contributions from eight birth metallicities and compared against the GWTC-5 constraints on the local rate (gray shaded region, which has been extended to $z=0.5$ for visualization). Three features of the cosmological BNS population can be inferred: (i) the local Universe ($z < 0.1$) is overwhelmingly dominated by progenitors born in the $0.1 - 2\,Z_\odot$ range; (ii) the contribution of metal-poor progenitors ($Z < 0.1\,Z_\odot$) peak at higher redshifts ($z\gtrsim3$); and finally, (iii) due to relatively short delay times, the the merger rate evolution broadly follows the cosmic star formation history, making BNS mergers potential tracers of cosmic chemical evolution.
Distributions of the total mass ($M_{\rm tot}$, left), chirp mass ($\mathcal{M}_{\rm c}$, center), and mass ratio ($q$, right) in \texttt{MODEL03} (top row) and \texttt{MODEL11} (bottom row), respectively. Both models use the same SN prescription \citep{Sukhbold+2016} but differ in their CE treatment (Table~\ref{tab:model_populations}). We show the cumulative populations merging within redshift intervals $z < 0.1$ (black), $z < 1$ (blue), and $z < 3$ (red), where each interval includes the previous ones. The vertical dashed lines and gray shaded regions represent the median and 90\% credible intervals for the events GW170817 (lower values) and GW190425 (higher values).
Caption Distributions of the total mass ($M_{\rm tot}$, left), chirp mass ($\mathcal{M}_{\rm c}$, center), and mass ratio ($q$, right) in \texttt{MODEL03} (top row) and \texttt{MODEL11} (bottom row), respectively. Both models use the same SN prescription \citep{Sukhbold+2016} but differ in their CE treatment (Table~\ref{tab:model_populations}). We show the cumulative populations merging within redshift intervals $z < 0.1$ (black), $z < 1$ (blue), and $z < 3$ (red), where each interval includes the previous ones. The vertical dashed lines and gray shaded regions represent the median and 90\% credible intervals for the events GW170817 (lower values) and GW190425 (higher values).
Distributions of the total mass ($M_{\rm tot}$, left), chirp mass ($\mathcal{M}_{\rm c}$, center), and mass ratio ($q$, right) in \texttt{MODEL03} (top row) and \texttt{MODEL11} (bottom row), respectively. Both models use the same SN prescription \citep{Sukhbold+2016} but differ in their CE treatment (Table~\ref{tab:model_populations}). We show the cumulative populations merging within redshift intervals $z < 0.1$ (black), $z < 1$ (blue), and $z < 3$ (red), where each interval includes the previous ones. The vertical dashed lines and gray shaded regions represent the median and 90\% credible intervals for the events GW170817 (lower values) and GW190425 (higher values).
Caption Distributions of the total mass ($M_{\rm tot}$, left), chirp mass ($\mathcal{M}_{\rm c}$, center), and mass ratio ($q$, right) in \texttt{MODEL03} (top row) and \texttt{MODEL11} (bottom row), respectively. Both models use the same SN prescription \citep{Sukhbold+2016} but differ in their CE treatment (Table~\ref{tab:model_populations}). We show the cumulative populations merging within redshift intervals $z < 0.1$ (black), $z < 1$ (blue), and $z < 3$ (red), where each interval includes the previous ones. The vertical dashed lines and gray shaded regions represent the median and 90\% credible intervals for the events GW170817 (lower values) and GW190425 (higher values).
Normalized probability densities of the pre-SN C/O core masses of the first-born (SN1) and second-born (SN2) NSs, for \texttt{MODEL03} (top two rows, using \citealt{Sukhbold+2016}) and \texttt{MODEL00} (bottom two rows, using \citealt{Patton&Sukhbold2020}), at $10^{-4}\,Z_\odot$ (left) and $\Zsun$ (right). The dashed histograms show all DNS systems formed, while solid histograms denote the merging subset. For both prescriptions, the merging systems at $10^{-4}\,Z_\odot$ are representative of the full population. At $Z_\odot$, only systems with more massive C/O cores ($\gtrsim 2.5\,M_\odot$) evolve into merging BNSs, revealing a metallicity-dependent selection effect. The pre-SN2 C/O core masses fall between $\sim 1.4$--$2\,M_\odot$ in all cases.
Caption Normalized probability densities of the pre-SN C/O core masses of the first-born (SN1) and second-born (SN2) NSs, for \texttt{MODEL03} (top two rows, using \citealt{Sukhbold+2016}) and \texttt{MODEL00} (bottom two rows, using \citealt{Patton&Sukhbold2020}), at $10^{-4}\,Z_\odot$ (left) and $\Zsun$ (right). The dashed histograms show all DNS systems formed, while solid histograms denote the merging subset. For both prescriptions, the merging systems at $10^{-4}\,Z_\odot$ are representative of the full population. At $Z_\odot$, only systems with more massive C/O cores ($\gtrsim 2.5\,M_\odot$) evolve into merging BNSs, revealing a metallicity-dependent selection effect. The pre-SN2 C/O core masses fall between $\sim 1.4$--$2\,M_\odot$ in all cases.
Normalized probability densities of the pre-SN C/O core masses of the first-born (SN1) and second-born (SN2) NSs, for \texttt{MODEL03} (top two rows, using \citealt{Sukhbold+2016}) and \texttt{MODEL00} (bottom two rows, using \citealt{Patton&Sukhbold2020}), at $10^{-4}\,Z_\odot$ (left) and $\Zsun$ (right). The dashed histograms show all DNS systems formed, while solid histograms denote the merging subset. For both prescriptions, the merging systems at $10^{-4}\,Z_\odot$ are representative of the full population. At $Z_\odot$, only systems with more massive C/O cores ($\gtrsim 2.5\,M_\odot$) evolve into merging BNSs, revealing a metallicity-dependent selection effect. The pre-SN2 C/O core masses fall between $\sim 1.4$--$2\,M_\odot$ in all cases.
Caption Normalized probability densities of the pre-SN C/O core masses of the first-born (SN1) and second-born (SN2) NSs, for \texttt{MODEL03} (top two rows, using \citealt{Sukhbold+2016}) and \texttt{MODEL00} (bottom two rows, using \citealt{Patton&Sukhbold2020}), at $10^{-4}\,Z_\odot$ (left) and $\Zsun$ (right). The dashed histograms show all DNS systems formed, while solid histograms denote the merging subset. For both prescriptions, the merging systems at $10^{-4}\,Z_\odot$ are representative of the full population. At $Z_\odot$, only systems with more massive C/O cores ($\gtrsim 2.5\,M_\odot$) evolve into merging BNSs, revealing a metallicity-dependent selection effect. The pre-SN2 C/O core masses fall between $\sim 1.4$--$2\,M_\odot$ in all cases.
Normalized mass ratio ($q$) distributions of merging BNSs in \texttt{MODEL03} at $10^{-4}\,Z_\odot$ (black) and $Z_\odot$ (gray). The $10^{-4}\,Z_\odot$ population produces more symmetric BNS mergers, with $\gtrsim 80\%$ having $q \ge 0.9$, compared to $\simeq 50\%$ at $\Zsun$, which also lacks the pronounced peak near equal masses at ($q \simeq 1$). While assessing precise NS mass values stretches the limits of standard analytical SN prescriptions, this systematic shift toward asymmetry at high metallicities is a robust, metallicity-driven trend arising from binary evolution.
Caption Normalized mass ratio ($q$) distributions of merging BNSs in \texttt{MODEL03} at $10^{-4}\,Z_\odot$ (black) and $Z_\odot$ (gray). The $10^{-4}\,Z_\odot$ population produces more symmetric BNS mergers, with $\gtrsim 80\%$ having $q \ge 0.9$, compared to $\simeq 50\%$ at $\Zsun$, which also lacks the pronounced peak near equal masses at ($q \simeq 1$). While assessing precise NS mass values stretches the limits of standard analytical SN prescriptions, this systematic shift toward asymmetry at high metallicities is a robust, metallicity-driven trend arising from binary evolution.
The cumulative number of detectable BNS mergers per observer year as a function of redshift for proposed 3G detectors, Cosmic Explorer (left) and the Einstein Telescope (right), is shown for all the population models considered in this work. The observationally calibrated models are highlighted in color and the remaining are shown in gray. (\textit{Top panels}) For Cosmic Explorer, solid and dashed lines show the 40\,km\,+\,20\,km and 20\,km\,+\,20\,km networks; for the Einstein Telescope, they show the 2L configuration with 15\,km arms and the triangle with 10\,km arms. (\textit{Bottom panels}) Ratio of the cumulative detections of the more sensitive configuration to that of the less sensitive one. Our observationally calibrated models predict $\simeq 10^4$--$10^5$ detectable BNS mergers per year out to a redshift $z\sim2$ (see Table \ref{tab:model_populations}).
Caption The cumulative number of detectable BNS mergers per observer year as a function of redshift for proposed 3G detectors, Cosmic Explorer (left) and the Einstein Telescope (right), is shown for all the population models considered in this work. The observationally calibrated models are highlighted in color and the remaining are shown in gray. (\textit{Top panels}) For Cosmic Explorer, solid and dashed lines show the 40\,km\,+\,20\,km and 20\,km\,+\,20\,km networks; for the Einstein Telescope, they show the 2L configuration with 15\,km arms and the triangle with 10\,km arms. (\textit{Bottom panels}) Ratio of the cumulative detections of the more sensitive configuration to that of the less sensitive one. Our observationally calibrated models predict $\simeq 10^4$--$10^5$ detectable BNS mergers per year out to a redshift $z\sim2$ (see Table \ref{tab:model_populations}).
The cumulative number of detectable BNS mergers per observer year as a function of redshift for proposed 3G detectors, Cosmic Explorer (left) and the Einstein Telescope (right), is shown for all the population models considered in this work. The observationally calibrated models are highlighted in color and the remaining are shown in gray. (\textit{Top panels}) For Cosmic Explorer, solid and dashed lines show the 40\,km\,+\,20\,km and 20\,km\,+\,20\,km networks; for the Einstein Telescope, they show the 2L configuration with 15\,km arms and the triangle with 10\,km arms. (\textit{Bottom panels}) Ratio of the cumulative detections of the more sensitive configuration to that of the less sensitive one. Our observationally calibrated models predict $\simeq 10^4$--$10^5$ detectable BNS mergers per year out to a redshift $z\sim2$ (see Table \ref{tab:model_populations}).
Caption The cumulative number of detectable BNS mergers per observer year as a function of redshift for proposed 3G detectors, Cosmic Explorer (left) and the Einstein Telescope (right), is shown for all the population models considered in this work. The observationally calibrated models are highlighted in color and the remaining are shown in gray. (\textit{Top panels}) For Cosmic Explorer, solid and dashed lines show the 40\,km\,+\,20\,km and 20\,km\,+\,20\,km networks; for the Einstein Telescope, they show the 2L configuration with 15\,km arms and the triangle with 10\,km arms. (\textit{Bottom panels}) Ratio of the cumulative detections of the more sensitive configuration to that of the less sensitive one. Our observationally calibrated models predict $\simeq 10^4$--$10^5$ detectable BNS mergers per year out to a redshift $z\sim2$ (see Table \ref{tab:model_populations}).
The total stellar mass formed per unit comoving volume integrated over cosmic history as a function of birth metallicity, calculated using Equation \ref{eq:sf_per_met}. We compare predictions from the cosmological simulation \texttt{IllustrisTNG} with three empirical prescriptions from \citet{Madau&Fragos2017}, \citet{Madau&Dickinson2014}\footnote{The \citet{Madau&Dickinson2014} star formation rate density, derived for a \citet{Salpeter1955} IMF, is multiplied by 0.66 to convert it to the \citet{Kroupa2001} IMF assumed in our populations \citep[following][]{Briel+2022}.} and \citet{Neijssel+2019} (only the first from this list and \texttt{IllustrisTNG} were adopted in this work). While all models agree that the bulk of stellar mass in the Universe is assembled at $Z > 10^{-2}\,Z_\odot$, they diverge drastically at lower metallicities. Specifically, \texttt{IllustrisTNG} predicts orders of magnitude more star formation at $Z \lesssim 10^{-2}\,Z_\odot$ compared to the sharp cutoffs inherent to the analytical models. The divergence at these low metallicity regimes constitutes a major systematic uncertainty when calculating the cosmological rates of BNS mergers originating from the most metal-poor progenitors.
Caption The total stellar mass formed per unit comoving volume integrated over cosmic history as a function of birth metallicity, calculated using Equation \ref{eq:sf_per_met}. We compare predictions from the cosmological simulation \texttt{IllustrisTNG} with three empirical prescriptions from \citet{Madau&Fragos2017}, \citet{Madau&Dickinson2014}\footnote{The \citet{Madau&Dickinson2014} star formation rate density, derived for a \citet{Salpeter1955} IMF, is multiplied by 0.66 to convert it to the \citet{Kroupa2001} IMF assumed in our populations \citep[following][]{Briel+2022}.} and \citet{Neijssel+2019} (only the first from this list and \texttt{IllustrisTNG} were adopted in this work). While all models agree that the bulk of stellar mass in the Universe is assembled at $Z > 10^{-2}\,Z_\odot$, they diverge drastically at lower metallicities. Specifically, \texttt{IllustrisTNG} predicts orders of magnitude more star formation at $Z \lesssim 10^{-2}\,Z_\odot$ compared to the sharp cutoffs inherent to the analytical models. The divergence at these low metallicity regimes constitutes a major systematic uncertainty when calculating the cosmological rates of BNS mergers originating from the most metal-poor progenitors.
References
  • Abac, A., et al. 2026, JCAP, 03, 081, doi: 10.1088/1475-7516/2026/03/081
  • Abbott, B. P., et al. 2017a, Class. Quant. Grav., 34, 044001, doi: 10.1088/1361-6382/aa51f4
  • Abbott, B. P., et al. 2017b, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  • Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Living Reviews in Relativity, 21, 3, doi: 10.1007/s41114-018-0012-9
  • Abbott, B. P., et al. 2020, Astrophys. J. Lett., 892, L3, doi: 10.3847/2041-8213/ab75f5
  • Aguilera-Dena, D. R., Langer, N., Antoniadis, J., et al.
  • 2022, A&A, 661, A60, doi: 10.1051/0004-6361/202142895
  • Aguilera-Dena, D. R., Müller, B., Antoniadis, J., et al.
  • 2023, A&A, 671, A134, doi: 10.1051/0004-6361/202243519
  • Andrews, J. J. 2020, ApJL, 900, L41, doi: 10.3847/2041-8213/abb1bf
  • Andrews, J. J., & Mandel, I. 2019, ApJL, 880, L8, doi: 10.3847/2041-8213/ab2ed1
  • Andrews, J. J., Bavera, S. S., Briel, M., et al. 2025, ApJS, 281, 3, doi: 10.3847/1538-4365/adfb78
  • 10−4
  • 10−3
  • 10−2
  • 10−1
  • 100
  • Metallicity Z/Z
  • 103
  • 104
  • 105
  • 106
  • 107
  • 108
  • Total stellar mass formed per unit comoving volume
  • [M
  • Mpc
  • −3
  • ]
  • IllustrisTNG
  • Neijssel+19
  • Madau+Fragos17
  • Madau+Dickinson14