The diffuse supernova neutrino background: a modern approach

Supernova neutrinos are a missing piece in the rapidly developing field of multimessenger astronomy. While the detection of a neutrino burst from an individual collapsing star remains a once-in-a-lifetime event — which occurred only once in recent history, with the historic observation of SN1987A Hirata (1987); Bionta et al. (1987); Alekseev et al. (1987) — there is confidence that, within a few years, the detection of supernova neutrinos will become commonplace, after the Diffuse Supernova Neutrino Background (DSNB) Bisnovatyi-Kogan and Seidov (1984); Krauss et al. (1984) starts to be observed at kiloton detectors Bays et al. (2012); Zhang et al. (2015); Abe et al. (2021); Abusleme et al. (2022); Harada et al. (2023); Cuesta Soria (2024).

Besides accelerating experimental progress, the DSNB is of great theoretical significance (see Beacom (2010); Lunardini (2016); Mirizzi et al. (2016); Ando et al. (2023); Tamborra (2025) for reviews), because it gives a global picture of the supernova population in the present universe and at cosmological distances as well. It is also an important probe of particle physics beyond the Standard Model, due to its unique sensitivity to new phenomena (e.g., neutrino decay, or neutrino-dark matter scattering) that develop over cosmological propagation distances (see Jeong et al. (2018); Das et al. (2024); Ivanez-Ballesteros and Volpe (2023); Balantekin et al. (2023); MacDonald et al. (2025) for recent examples).

Predicting the DSNB is challenging due to the myriad of effects that influence it; the main ones being stellar evolution, the physics of collapse simulations, the core-collapse rate history, and neutrino physics. Each element needs to be modeled in detail across the population of progenitor stars, which is currently unfeasible.

Early computations of the DSNB  – e.g., Woosley et al. (1986); Malaney (1997); Hartmann and Woosley (1997); Lunardini (2006); Horiuchi et al. (2009) – were simplified, using a single parameterized supernova neutrino spectrum, combined with the cosmological star formation history. Later predictions developed more detail, incorporating the variation of neutrino emission across the population in the form of dependence on the the zero-age main sequence (ZAMS) mass (Totani and Sato, 1995; Lunardini and Tamborra, 2012; Horiuchi et al., 2018; Priya and Lunardini, 2017; Kresse et al., 2021; Ekanger et al., 2024) and including black-hole-forming collapses (see, e.g., (Lunardini, 2009; Mathews et al., 2014; Nakazato et al., 2015; Ashida and Nakazato, 2022; Ashida et al., 2023; Nakazato et al., 2024; Martínez-Miravé et al., 2024)). These works used the results of 1D and 2D collapse simulations (e.g, Fischer et al., 2010; Hüdepohl, 2013; Nakazato, 2013; Mirizzi et al., 2016; Vartanyan and Burrows, 2023), which were obtained for certain sets of progenitors where stellar evolution effects were only partially included. Binary evolution effects were recognized as potentially important ¹¹1Pioneering works of Vartanyan et al. (2021); Wang and Pan (2024); Schneider et al. (2021) employ some binary effects in the modeling of a supernova neutrino burst. The binary progenitor models are computed in Ref. (Laplace et al., 2021; Schneider et al., 2024). Some of the data is publically avialble in https://doi.org/10.5281/zenodo.14959596. — since most massive stars are in binary systems Sana et al. (2012) — but were not included in DSNB calculations. An exception is the work by Kresse, Ertl and Janka Kresse et al. (2021), where binary-evolved stars were schematically represented by a set of helium stars Woosley (2019).

Over time, numerical mappings of the core collapse parameter space have led to identifying some key parameters - other than the ZAMS mass - that most influence the neutrino emission and connect it with stellar evolution: the mass of the Carbon-Oxygen (CO) core and the compactness parameter Limongi and Chieffi (2018); Takahashi et al. (2023). In particular, there is now consensus that the neutrino emission depends mainly on compactness O’Connor and Ott (2013); Nakamura et al. (2015); Warren et al. (2020), and the compactness is roughly determined primarily by the CO core mass at the advanced stage (see Takahashi et al., 2023; Patton and Sukhbold, 2020, for the detail), which is a typical outcome of stellar evolution simulations. The immediate implication of these advances is the possibility to predict the neutrino emission for a given model of stellar population if its CO core mass distribution is known Horiuchi et al. (2018). This opens up the possibility to consistently interface two important problems: the modeling of neutrino emission from collapsing stars, and the way supernova progenitors evolve over their lifetime including binary effects.

This idea was first applied to the DSNB in the exploratory studies by Horiuchi et al. (2021) (Horiuchi et al., 2021, hereafter H21), with a focus on the effect of binary evolution. That work is extension of Horiuchi et al. (2018) (Horiuchi et al., 2018, hereafter H18), where over 100 progenitors from Ref. Nakamura et al. (2015) was used in determining the DSNB but without binary considerations. In H21, 20 progenitors were used from the study of Ref. Summa et al. (2016) and binary effects were considered. Note that most, though not all, of the core-collapse simulations employed in H21 typically terminate in a short time, $t<1\,{\rm s}$ Nakamura et al. (2015); Summa et al. (2016), and their neutrino flux outputs are extrapolated to obtain the whole time evolution of the neutrino emission Arcones et al. (2007). This introduces a systematic uncertainty Ekanger et al. (2022) which can be updated Ekanger et al. (2024) using long-term simulations. In H21, binary effects were modeled using the population synthesis method developed by Ref. Hurley et al. (2002) and updated in Ref. Kinugawa et al. (2014) (described in more detail below). The modifications to the CO mass distribution due to binary interactions were quantified using six binary population synthesis realizations, but since the number of supernova progenitors used to study neutrino emission was sparse (20 progenitors), a schematic mapping was developed: the neutrino emission was modeled using the core compactness parameter, and the mapping between core compactness and CO core mass was quantified. This allowed H21 to model the neutrino emission from hundreds of thousands of supernova progenitors in the binary population synthesis models. See Fig. 1 for an illustration of the procedure.

Here we further develop the method of H21. We discuss its elements in detail, and to what degree they influence the DSNB. For each element, we employ the most recent, state-of-the-art predictions. Among these, an important update is the neutrino emission parameters (spectra and luminosities) from the recent set of 100 multi-second 2D core collapse simulations of Vartanyan and Burrows (2023, hereafter VB23). We also use the results of two stellar evolution codes, BPASS Patton et al. (2022) and the code adopted in H21 by Kinugawa, Inayoshi, Hotokezaka, Nakauchi, and Nakamura (Kinugawa et al., 2014, KIHNN from here on), where binary evolution effects are modeled in detail, including key processes such as Roche lobe overflow, common envelope evolution, and stellar mergers. These effects can significantly change the final fate of massive stars: for instance, mass transfer can change the stellar masses, and can strip hydrogen envelopes, leading to helium stars or Wolf-Rayet-like progenitors, while stellar mergers can produce more massive cores than single-star evolution would allow. Such processes alter the CO core mass and compactness at collapse, thereby influencing the resulting neutrino emission. The BPASS code and the KIHNN code track these effects across a wide range of initial binary parameters and assumptions of binary interaction parameters. Our study is meant to lay out a methodology that is particularly well suited for future developments. Indeed, once substantial new sets of numerical simulations are accessible, it will be feasible to derive the fundamental dependencies on CO-core mass and compactness from them, and subsequently update the DSNB calculation in a consistent manner.

The paper is structured as follows. In Sec. II, the method is introduced, along with some basic general information. In Sec. III, the neutrino emission from collapsing stars is described in detail, with emphasis on its dependence on compactness. Sec. IV contains a discussion of binary evolution effects. In Sec. V, the elements presented in the previous sections are combined into the calculation of the DSNB, and results are presented. A summary and discussion follows in Sec. VI.

Before delving into a detailed description, here we summarize the main parts of this work and clarify their interplay. We also define the key quantities and notation.

To calculate the DSNB, we require: (i) the neutrino flux of individual core-collapses, and (ii) the distribution of the population of these events. The neutrino flux is given by core-collapse simulations with neutrino radiation transfer and the distribution of the event is provided by stellar evolution calculations. To bridge the two sets of simulations, we have to parametrize the simulations results using compactness and CO-core mass. We illustrate the procedure in Fig. 1. The elements shown in the figure are briefly described below.

To model the neutrino emission from a supernova, we employ the results of the most recent, state-of-the-art models: 100 core-collapse simulations developed by VB23, which are two-dimensional (2D), and extend up to 4 s post-bounce. In these simulations, black hole formation appears for ZAMS mass $M\simeq$ 12–15 $\mathrel{{M_{\odot}}}$, and this trend is different from the previously established picture, where black holes were formed for $M\geq 20\,\mathrel{{M_{\odot}}}$ (e.g., Woosley et al., 2002). To characterize the neutrino flux, we use the compactness parameter O’Connor and Ott (2011), defined as:

which is evaluated at the pre-collapse phase. This quantity is obtained by integrating the mass density up to a certain radius, $R$, to obtain the mass enclosed within that radius, $M_{\rm e}$. It is customary to fix the value of $M_{\rm e}$; here we adopt $M_{\rm e}=2.5\mathrel{{M_{\odot}}}$, following H21. Note that other choices of $M_{e}$ are correlated with that of $M_{\rm e}=2.5\mathrel{{M_{\odot}}}$ Pejcha and Thompson (2015), and lead to very similar results in our study. The compactness parameter has the merit of being simple, and it is widely used to predict the key observables, such as compact remnant mass, explosion energy and the properties of the neutrino emission Takahashi et al. (2023); O’Connor and Ott (2013); Nakamura et al. (2015); Warren et al. (2020).²²2However, we note that compactness alone may not be the best predictor of supernova explodability (e.g., Wang et al., 2022; Tsang et al., 2022; Burrows et al., 2024). We refer to Sec. III.2 for details.

To portray the distributions of the core-collapse events, the mass of the CO core, $\mathrel{{M_{\rm CO}}}$, is a good parameter of choice. The reason is practicality: in large stellar evolution grids, evolving all models to the pre-collapse phase is computationally challenging, and obtaining the compactness is difficult. Instead, $\mathrel{{M_{\rm CO}}}$ is more readily available in the most stellar evolution studies (e.g., Patton et al., 2022; Horiuchi et al., 2021; Fragos et al., 2023). This motivates using a mapping to connect $\mathrel{{\xi_{2.5}}}$ and $\mathrel{{M_{\rm CO}}}$ to connect stellar evolution simulations and core collapse ones (Fig. 1; see Sec. III.1 for details). Rigorously, the value of $\mathrel{{M_{\rm CO}}}$ at earlier stages, typically C or Ne ignition — which is most often provided by stellar evolution simulations — may slightly differ from that at the onset of core collapse.³³3In BPASS, stars are evolved until core-collapse, so this issue does not apply. Also, the compactness depends on the composition of the CO-core Patton and Sukhbold (2020) and this effect is not included in our work. These problems should be addressed in future studies.

A novel aspect of our study lies in the detailed comparison between single and binary star evolution, according to different numerical models. In Sec. IV, we present the distributions of $\mathrel{{M_{\rm CO}}}$ based on two recent stellar evolution models: BPASS Patton et al. (2022) and the KIHNN model Kinugawa et al. (2014), which was employed in H21. These distributions are obtained by assuming certain distribution of initial binary parameters, such as mass ratio, orbital period (separation), and eccentricity.

Considering the cosmological time evolution and effect of neutrino oscillation, we obtain DSNB flux, as described in Sec. V.

Binary star evolution can involve a wide range of physical processes that are absent in single star evolution, and these processes can significantly alter the outcome of stellar evolution. One of the most common is mass transfer, in which material flows from one star to its companion, often through Roche lobe overflow. Another key process is common envelope evolution, where both stars share a single, extended envelope and spiral toward each other, potentially leading to a merger. Stellar mergers can produce unusually massive stars with higher core masses than possible via isolated evolution. Binary interactions can also affect angular momentum, rotation, and core structure, which in turn influence the final compact remnant and its neutrino emission.

Because binary evolution proceeds differently depending on the initial parameters of the system, such as the primary mass, mass ratio, orbital separation, and eccentricity, it is necessary to investigate these outcomes statistically by sampling from distributions of these initial parameters. This approach is known as binary population synthesis, and it enables a systematic study of the diverse evolutionary paths and their effects on observable quantities such as the diffuse supernova neutrino background. In addition, binary interactions are influenced by physical processes whose outcomes remain uncertain, such as the accretion efficiency of mass transfer ($\beta$) and the efficiency of common envelope ejection ($\alpha\lambda$), making parameterized treatments necessary within population synthesis calculations.

We adopt two binary population synthesis results. One is the BPASS v2.2 model Patton et al. (2022) and the other is the KIHNN  model Kinugawa et al. (2014). In the BPASS runs, a total mass $M_{\rm tot}=10^{6}\mathrel{{M_{\odot}}}$ of single stars, and an equal mass of binary systems are considered. BPASS is a hybrid Binary Population Synthesis (BPS) code. It evolves its single stars and the primaries of binary stars with a detailed stellar evolution code, while approximating the evolution of the secondaries using the semi- analytic evolution equations from Hurley et al. (2002). It then evolves the secondaries in detail once the primary star’s evolution terminates.

A detailed description of BPASS is given in Stanway and Eldridge (2018). In the single star models, the ZAMS mass range span $0.2$–$100\mathrel{{M_{\odot}}}$ in steps of $0.1\mathrel{{M_{\odot}}}$; in the binary models masses span $4.5$–$100\mathrel{{M_{\odot}}}$ for the primary star’s mass, $M_{1}$, with increments of $0.5\mathrel{{M_{\odot}}}$. For the binary systems, the mass ratio, $q$, varies between 0.1 and 0.9 in increments of 0.1, and the orbital period, $P$, ranges between $10^{1.4}$ and $10^{4}$ days in steps of 0.2 dex in log space. The distributions of $q$ and $P$ follows Moe and Di Stefano (2017). All models are for solar metallicity. Stars are evolved until core carbon ignition. BPASS assumes a Kroupa (2001) initial mass function (IMF) for masses between 0.5 and $1\mathrel{{M_{\odot}}}$, and a Salpeter (1955) IMF (which follows $M^{-2.35}$) (Salpeter, 1955) for stars with masses heavier than $1\mathrel{{M_{\odot}}}$. The mass and radius of the CO-core are computed by fixing the boundary of the CO-core at the radius where the helium mass fraction equals 0.1.

The KIHNN model adopts the fiducial model ($\alpha\lambda=1$) from H21 and uses the Abt_MT1_CE1 model as described in Kinugawa et al. (2024).⁷⁷7The fiducial model ($\alpha\lambda=1$) in H21 assumes a primary star mass range of 3 $M_{\odot}$ to 140 $M_{\odot}$, while the Abt_MT1_CE1 model in Kinugawa et al. (2024) assumes 3 $M_{\odot}$ to 100 $M_{\odot}$. However, since the IMF is Salpeter, the contribution from stars in the 100 $M_{\odot}$ to 140 $M_{\odot}$ range is negligible, making the two models effectively equivalent. In this study, we adopt the range of 3 $M_{\odot}$ to 100 $M_{\odot}$. The KIHNN model is calculated by a population synthesis code based on an updated version of the binary stars evolution code by Hurley et al. (2002), with details described in Ref. Kinugawa et al. (2014); Horiuchi et al. (2021); Kinugawa et al. (2024). The main differences include the treatment of mass transfer stability, the handling of the common envelope phase, and the evolution of the cores of rotating stars.

In the KIHNN model, we assume solar metallicity and adopt a Salpeter IMF for the mass of primary stars in the range of $3$ to $100\,M_{\odot}$, a constant initial mass ratio function (IMRF) (Kobulnicky et al., 2007), a logarithmically uniform initial separation function (ISF) proportional to the inverse of the separation, $1/a$, (Abt, 1983), and an initial eccentricity function (IEF) proportional to the eccentricity, $e$ (Heggie, 1975). In addition, pulsar kick velocities are modeled using a Maxwellian distribution with a dispersion of $\sigma=250\,\mathrm{km\ s}^{-1}$ (Hobbs et al., 2005). In the KIHNN model, we use binary parameters set as $(\beta,\alpha\lambda)=(1,1)$. Here, $\beta=1$ signifies that no mass is lost from the binary system during the mass transfer. The quantity $\alpha\lambda$ are the parameters of Common Envelope (CE) phase, where $\alpha$ is the efficiency with which orbital energy is used to unbind the envelope, and $\lambda$ characterizes the binding energy of the envelope. In the KIHNN model, $\alpha\lambda$ is treated as a single parameter. This parameterization plays a crucial role in determining whether the binary survives the common envelope phase or merges prematurely; see Horiuchi et al. (2021) for details. When the accretor is a compact object, accretion is limited by the Eddington limit, and any excess mass is ejected from the system. If a stellar merger occurs due to a common envelope phase before the star becomes a compact object, the resulting star is treated as a highly rotating star. In this case, the increase in core mass due to rotation is taken into account Horiuchi et al. (2021); Kinugawa et al. (2024). A total of $10^{5}$ binary systems were simulated using this configuration.

Fig. 4 and 5 show the CO core distributions obtained by the two models, for single and binary stellar evolution. Despite the differences that exist between different simulations, some predictions are robust. One of them is that, as a result of binary case, stars are produced with high $\rm CO$ core mass, $\mathrel{{M_{\rm CO}}}\gtrsim 10\mathrel{{M_{\odot}}}$, that could not be realized for single star evolution.⁸⁸8Note that this feature depends on the treatment of mass loss, see Figure 4 of Fragos et al. (2023) and Figure 1 of Patton et al. (2022) These fall in the region of medium-high compactness (see Fig. 2 (a)), and therefore have a larger neutrino emission (high $E_{\rm tot}$) than most neutron-star-forming single-evolved stars. In particular, binary-evolved stars populate the (possible) second peak of the $\mathrel{{\xi_{2.5}}}$ distribution ($\mathrel{{M_{\rm CO}}}\sim 12\mathrel{{M_{\odot}}}$, see Fig. 2 (b)), where $\mathrel{{\bar{\nu}}_{e}}$s could be emitted with a factor of $\sim 2$ larger total energy ($E_{\rm tot,\mathrel{{\bar{\nu}}_{e}}}\simeq 1.2\times 10^{53}\,{\rm erg}$) and a $\sim 10\%$ higher average energy than for a NSFC  with $\mathrel{{M_{\rm CO}}}\simeq 2.5\mathrel{{M_{\odot}}}$. These progenitors with $\mathrel{{M_{\rm CO}}}\gtrsim 10\mathrel{{M_{\odot}}}$ are exceedingly rare, however, amounting to $\sim 1\%$ or less of the entire population.

Another prediction is a wide distribution in $\mathrel{{M_{\rm CO}}}$ that becomes possible from binary evolution even for progenitors of relatively modest ZAMS mass. For example, a binary system where $M_{1}<M_{2}\simeq 15\mathrel{{M_{\odot}}}$ – for which single evolution predicts $\mathrel{{M_{\rm CO}}}\lesssim 4\mathrel{{M_{\odot}}}$– could lead to the formation of a CO core with up to $\mathrel{{M_{\rm CO}}}\simeq 10\mathrel{{M_{\odot}}}$. Interestingly, however, the distributions of the stellar population with $\mathrel{{M_{\rm CO}}}$ are relatively similar for the two cases of single and binary evolutions. For the BPASS results, the main difference is that for binary evolution the local maximum at $\mathrel{{M_{\rm CO}}}\sim 8$–$10\mathrel{{M_{\odot}}}$ becomes de-populated, whereas the distribution remains practically the same as the single evolution case in the region $\mathrel{{M_{\rm CO}}}\lesssim 4\mathrel{{M_{\odot}}}$. In the KIHNN model, instead, a slight overpopulation of the low $\mathrel{{M_{\rm CO}}}$ region is found for binary evolution.

We now combine the results of the previous sections with the cosmological rate of core collapse, to compute the DSNB.

For a given population synthesis model (BPASS or KIHNN model), the first step is to compute the following fraction:

where the index $i$ indicates the case of single or binary star evolution ($i=S,B$) and $F_{w}(E,\mathrel{{M_{\rm CO}}})$ is the individual star neutrino emission in the species $w$, which is derived by combining Eq. (2) with the function $\mathrel{{\xi_{2.5}}}=\mathrel{{\xi_{2.5}}}\!(\mathrel{{M_{\rm CO}}})$ as shown in Fig. 2 (b). The quantity $\left(\frac{dN}{d\mathrel{{M_{\rm CO}}}}\right)_{i}$ is the distribution of stars in $\mathrel{{M_{\rm CO}}}$, which is displayed in Fig. 4 (c), and Fig. 5 (c) for BPASS and KIHNN model, respctively.

For $i=S$, the expression in Eq. (3) is the population-averaged neutrino flux, whereas for $i=B$, it is the average multiplied by a normalization factor that accounts for binary evolution effects changing the total number of stars that undergo core collapse. We study two cases: the first is one where all the collapses form neutron stars, therefore $\mathrel{{\xi_{2.5}}}=\mathrel{{\xi_{2.5}}}\!(\mathrel{{M_{\rm CO}}})$ as shown in the curves in Fig. 3. The second case is with black hole formation in the range $2.0<\mathrel{{M_{\rm CO}}}/M_{\odot}<3.1$ (see Sec. III.2), which corresponds to $\sim 20\%$ of all collapses forming black holes. Here, the neutrino spectrum parameters are taken to be fixed to their BHFC values (empty markers in Fig. 3) in the interval $2.0<\mathrel{{M_{\rm CO}}}/M_{\odot}<3.1$, and to the NSFC curves in the complementary region of $\mathrel{{M_{\rm CO}}}$.

Following H21, we make the assumption that $\left(\frac{dN}{d\mathrel{{M_{\rm CO}}}}\right)_{i}$ — and therefore $f^{w}_{i}(E)$ — does not depend on the redshift $z$. This is justified by the fact that binary evolution effects occur on a time scale which is much shorter than the lifetime of a supernova progenitor, even when binary effects are included.

For the cosmological rate of core collapses, we use the parametrization in Yuksel et al. (2008) (which was also used in H21):

where $\dot{\rho}_{0}=1.3\times 10^{-4}\,\text{yr}^{-1}\,\text{Mpc}^{-3}$, $a=3.4$, $b=-0.3$, $c=-3.5$, $B=5000$, $C=9$ and $\eta=-10$. We checked that using other functional forms (e.g., the ones in Ekanger et al. (2024)) leads to nearly identical results if the normalization $\dot{\rho}_{0}$ is the same.

From the quantities above, one obtains the DSNB for a single neutrino species without flavor oscillations:

where $E^{\prime}=E(1+z)$; the fractions of cosmological energy density in matter and dark energy are $\Omega_{\rm m}=0.3$ and $\Omega_{\Lambda}=0.7$, and the Hubble constant is fixed to $H_{0}=71\,{\rm km}\,{\rm s}^{-1}{\rm Mpc}^{-1}$. We fixed $z_{\rm max}=4.5$; due to the fast decline of the supernova rate at $z\gtrsim 4$, the result of the integral depends only minimally on the value of this parameter. Here $c$ in the equation is the speed of light.

To compute the realistic flux of $\mathrel{{\bar{\nu}}_{e}}$ and $\mathrel{{\nu_{e}}}$ in a detector on Earth, one must include the effect of neutrino flavor oscillation. This is done by introducing an effective $\mathrel{{\bar{\nu}}_{e}}$ survival probability, ${\bar{p}}$, and writing the (oscillated) diffuse $\mathrel{{\bar{\nu}}_{e}}$ flux as:

A similar expression as Eq. (6) holds for $\mathrel{{\nu_{e}}}$, with the replacement $\bar{p}\rightarrow p$. The parameters $\bar{p}$ and $p$ can be computed by modeling the matter-driven and neutrino-driven (i.e., collective flavor oscillations) inside the star, as well as oscillations inside the Earth. Their values are uncertain, due to the uncertainty in the treatment of collective flavor oscillations. For our purpose, it is sufficient to approximate them as constants, with values falling in intervals $\bar{p}=0$–$0.7$ and $p=0$–$0.3$. The extremes of these intervals are (approximately) the predictions of matter-driven conversion depending on the choice of neutrino mass ordering (specifically, $\bar{p}\simeq\sin^{2}\theta_{13}$–$\cos^{2}\theta_{12}$ and $p\simeq\sin^{2}\theta_{13}$–$\sin^{2}\theta_{12}$ see, e.g., Lunardini and Tamborra (2012) for a discussion and further references).

Fig. 6 illustrates the effect of including binary evolution effects and/or black hole formation on the DSNB, for the un-oscillated flavor fluxes and the two simulations of binary evolution.

It appears that binary evolution effects are moderate: for the KIHNN model, they amount to a $\sim 10-15\%$ increase in the normalization of the flux, reflecting the increase in the number of collapse candidates — especially at low $\mathrel{{M_{\rm CO}}}$ — compared to single-star evolution. That is consistent with the results of H21, and different from Kresse et al. (2021), who found that including binary effects decreases the DSNB, since mass stripping decreases the CO core mass Woosley (2019). Mergers of low-mass stars, whose individual masses are below the core-collapse threshold, can increase the DSNB if the combined mass exceeds the threshold. In contrast, mergers between stars already massive enough to undergo core collapse individually tend to reduce the DSNB. The net effect of binary interactions is thus determined by these competing processes (see the discussion in H21, ).

For the case of BPASS, binary evolution effect are at the level of $\sim-4\%$, and are invisible in Fig. 6. The difference in the DSNB results between the cases with BPASS and with the KIHNN model can partially be traced to the assumptions of the evolution of the core after a stellar merger. In the KIHNN model, the merged remnants are assumed to be highly rotating. The high rotation is then considered to lead to an increase in the core mass during the subsequent evolution.

The effect of including black hole formation leads to a high energy tail in the (un-oscillated) $\mathrel{{\bar{\nu}}_{e}}$ and $\mathrel{{\nu_{e}}}$ spectra, which is consistent with previous literature (e.g., Lunardini (2009)). The flux increase — with respect to NSFC only — is of $\sim 50\%$ at $E=35$ MeV. For $\mathrel{{\nu_{x}}}$ there is an increase of similar magnitude in the normalization, while the energy spectrum is practically unmodified, due to the very modest difference in $\langle E\rangle_{x}$ between NSFC and BHFC (Fig. 3).

In Fig. 7, the $\mathrel{{\bar{\nu}}_{e}}$ and $\mathrel{{\nu_{e}}}$ components of the DSNB, $\Phi^{\bar{e},i}(E)$ and $\Phi^{e,i}(E)$ (which include flavor oscillations), are shown, for different values of the flavor-swapping parameters. Binary evolution and black hole formation are included. For both species, the main effect of flavor oscillations is to decrease the expected flux, compared to the no-oscillation case. The decrease is largest for $\mathrel{{\nu_{e}}}$ near the peak of the flux, $E\sim 3$ MeV, where it can amount up to a factor of $\sim 2$; it is at the level of $\sim 20\%$ or so in the $\mathrel{{\bar{\nu}}_{e}}$ channel. This trend is in contrast with typical results that can be found in the literature, where oscillations tend to increase the $\mathrel{{\nu_{e}}}$ and $\mathrel{{\bar{\nu}}_{e}}$ fluxes, at least in the region of sensitivity of current detectors. The difference is due to the fact that here the originally produced $\mathrel{{\nu_{x}}}$ (part of which is converted into $\mathrel{{\nu_{e}}}$ and $\mathrel{{\bar{\nu}}_{e}}$) have lower luminosity than the other species, and spectra that are very similar to those for $\mathrel{{\bar{\nu}}_{e}}$. Instead, most literature assume $\mathrel{{\nu_{x}}}$ fluxes with (approximately) equal luminosity and harder spectra than the electron flavors.

We have developed a methodology for modeling the DSNB that is based on the mapping from CO core mass (at advanced pre-supernova stage) to compactness, as summarized in Figs. 1 and 2. The method is suitable for the systematic inclusion of (i) neutrino spectra from large sets of numerical core collapse simulations, and (ii) population synthesis models. We have applied it to state-of-the-art models; in particular, we have processed the results of the extended suite of simulations from Vartanyan and Burrows (2023), and we have identified the relevant dependence of the neutrino spectra and fluxes on the compactness parameter. We have also analyzed the results of two modern population synthesis codes, BPASS Patton et al. (2022) and the KIHNN model Kinugawa et al. (2014), and obtained the distribution of the stellar populations with respect to the CO-core mass, with and without binary evolution effects.

We can now draw a number of conclusions of our work:

1. 1.

   our analysis of the numerical results of Vartanyan and Burrows (2023) for neutron-star-forming collapses shows a dependence of the neutrino flux parameters on $\mathrel{{\xi_{2.5}}}$ which is linear for the total energies emitted (per flavor) and like $\mathrel{{\xi_{2.5}}}^{1/3}$ for the average energies and rms energies. The total energy emitted in each of the non-electron flavors is lower by $\sim 20$–$50\%$ than the electron flavors. For antineutrinos, this is actually the main difference between $\mathrel{{\bar{\nu}}_{e}}$ and $\mathrel{{\nu_{x}}}$, because the spectra of these two species are very similar. The spectrum of $\mathrel{{\nu_{e}}}$ is colder than that of other species (a $\sim 30\%$ difference in the average energies) which consistent with older models.

2. 2.

   The impact of neutrino flavor oscillations on the DSNB is smaller than in most previous literature, due to electron and non-electron flavors having similar spectra at the production site, especially in the antineutrino channel. In contrast with the traditional expectation, oscillations could mostly modify the flux normalization rather than the energy spectrum, in the energy window of sensitivity of current and near future experiments ($E\gtrsim 10$ MeV). In the neutrino channel, flavor conversion can result in harder $\mathrel{{\nu_{e}}}$ spectra, mainly in the form of a reduction of the flux peak at $E\sim 4$ MeV.

3. 3.

   In the work of Vartanyan and Burrows (2023), black-hole-forming collapses are predicted to occur in $\sim 20\%$ of all collapses, for $2.0<\mathrel{{M_{\rm CO}}}/M_{\odot}<3.1$ (meaning medium mass progenitors in case of single-evolving stars, $12\leq M/M_{\odot}\leq 15.38$). These black-hole forming collapses have more energetic neutrino spectra and higher luminosities, and therefore increase the DSNB in its high energy tail and in its overall normalization. The effect is of $\sim 40\%$–$60\%$ size.

4. 4.

   due to the modest dependence of the neutrino spectral parameters on $\mathrel{{\xi_{2.5}}}$ — and therefore on $\mathrel{{M_{\rm CO}}}$ — the main effect of binary evolution on the DSNB is a moderate change in its normalization. The change is largest, up to a $\sim 15\%$ increase, for the KIHNN model, and is caused mostly by the net increase in the number of collapsing stars. It is interesting that binary population synthesis predicts the existence of stars with more massive CO cores, $\mathrel{{M_{\rm CO}}}\gtrsim 10\mathrel{{M_{\odot}}}$, for which more energetic neutrino emission (than most single-evolved stars) is expected. However, these high $\mathrel{{M_{\rm CO}}}$ stars amount to only $\sim 1\%$ of the population; this fact, combined with weak dependence of the neutrino spectral parameters with $\mathrel{{M_{\rm CO}}}$ results in a negligible effect on the DSNB.

Before concluding, a number of cautionary remarks are in order. One of them is that our results suffer from uncertainties related to the scatter (see, e.g., Fig. 2 (b)) and sparseness (e.g., Fig. 5, central pane) of the numerical results that we have used. For this reason, we chose to consider small numerical features as not significant, and centered the discussion on then more robust, general trends that we found.

The path and rate of black hole formation remain highly uncertain, and these uncertainties could have a significant impact on the DSNB. In the core collapse simulations we used, a significant fraction of models exhibit long-lived proto neutron stars, $\geq 3.5\,{\rm s}$, without driving the supernova explosion. Typically, however, in the case that the explosion is not initiated, black hole formation occurs within $\mathcal{O}(1\,{\rm s})$ O’Connor and Ott (2011), with the timescale depending on the equation of state Liebendoerfer et al. (2004); Sumiyoshi et al. (2007); O’Connor and Ott (2011); Nakazato (2013); Walk et al. (2020). In some models, the black hole formation could occur within an exploding model Chan et al. (2018); Pan et al. (2018); Kuroda et al. (2018); Burrows et al. (2023). Mass accretion — enhancing the neutrino emission — could continue even after the shock revival, with the timescale of $\mathcal{O}(10\,{\rm s})$ and have large influence on DSNB Nakazato et al. (2024); Akaho et al. (2024). Future studies will require an extensive set of 3D simulations Bollig et al. (2021); Nakamura et al. (2024); Wang and Burrows (2023) to explore these formation scenarios in great detail.

Some uncertainties are also associated with neutrino flavor conversion. Recently, neutrino oscillations induced by the neutrino-neutrino interactions have been studied in detail Duan et al. (2010); Mirizzi et al. (2016); Tamborra and Shalgar (2020); Capozzi and Saviano (2022); Richers and Sen (2022); Volpe (2023). Though self-consistent simulations in several time snapshots have been performed employing this effects (Nagakura and Zaizen, 2023; Xiong et al., 2024), we have not found a simple way to incorporate this effect accurately in the global simulations (however see Ehring et al., 2023; Mori et al., 2025; Wang and Burrows, 2025, for the phenomenological modeling). We have to wait for the establishment of an appropriate treatment of collective oscillations to accurately assess their impact on the DSNB.

The study of the evolution of binary systems has yet to reach maturity Tauris and van den Heuvel (2023); Marchant and Bodensteiner (2024), and therefore current descriptions have significant uncertainties. It is possible that, as more detailed treatments are developed, a larger impact of binary evolution on the DSNB may emerge. For example, including possible highly-rotating stars may enhance the mixing of compositions and produce larger CO-core masses Horiuchi et al. (2021), that will result in higher neutrino emission. Perhaps new explosion channels incorporating rotation and magnetic fields can appear Summa et al. (2018); Takiwaki et al. (2021); Matsumoto et al. (2022, 2023); Obergaulinger and Aloy (2022); Shibagaki et al. (2024); Powell and Müller (2024); Kuroda and Shibata (2024). Also more sophisticated treatment of stellar mergers could alter the core mass distribution.

It is important to note that the effects we studied here are comparable or smaller than the uncertainty on the cosmological rate of core collapse, which is several tens of per cent or larger. Therefore, it may not possible to test them in the near future, when the first DSNB data become available. Indirectly, our study contributes to stress the importance of working to reduce the uncertainty on the core collapse rate.

In closing, we have produced a new, state-of-the-art model for the DSNB, which includes updated neutrino fluxes and spectra, and effects of binary evolution. We have done so using a modern method that lends itself to further developments, especially those that will be driven by large sets of numerical simulations. Our analysis highlights the potential of the DSNB as a test of a wide variety of physical phenomena — ranging from the microphysics inside collapsing stars, to the way stars that are born with companions evolve over their lifetime — and how such potential requires reducing the numerous uncertainties that are present in order to be fulfilled.