Add Diffuse Supernova Neutrino Background (DSNB) flux calculation module

[dristii-32]

Summary

This PR adds a new snewpy.dsnb module for computing the Diffuse Supernova Neutrino Background (DSNB) electron antineutrino flux and inverse beta decay (IBD) event rates.

The implementation follows:

• Li, Vagins & Wurm (2022), arXiv:2201.12920
• Hopkins & Beacom (2006) core-collapse supernova rate
• Keil, Raffelt & Janka (2003) pinched neutrino spectra
• Strumia & Vissani (2003) inverse beta decay cross section

Features

• Redshift-dependent core-collapse supernova rate
• Pinched supernova neutrino emission spectra
• Vectorized cosmological DSNB flux integration
• IBD event-rate calculation
• Optional detector energy-resolution smearing

Testing

The module imports successfully within SNEWPY.

Added a dedicated test suite:

pytest python/snewpy/test/test_dsnb.py

All tests pass:

27 passed

Notes

This contribution is based on an independently developed and tested implementation that has also been released as a standalone Python package on PyPI. The code has been adapted here for integration into SNEWPY.

[JostMigenda]

Hello Dristii, welcome to SNEWPY!

Before we go into the details of this code, my first question is at a much higher level: As far as I can tell, this code has really no connection to the existing SNEWPY package, right? (Apart from being in the same area of astrophysics, of course.) So to me, it's not at all clear what advantages there would be to including this as part of SNEWPY?

[dristii-32]

Thanks for the welcome and the question.

The motivation for this was that the DSNB is the natural diffuse equivalent of the time dependent supernova neutrino fluxes already supported by SNEWPY.

SNEWPY currently includes interfaces to supernova neutrino emission models and tools for computing the observable fluxes and event rates from individual core-collapse supernovae. But many neutrino experiments (Super-Kamiokande with Gd, JUNO, future detectors etc) also look at the diffuse supernova neutrino background, which is obtained by summing the contribution from the whole cosmic population of core-collapse supernovae.

The goal of this module is thus to expand SNEWPY from single-supernova calculations to the cosmological DSNB, so that users can stay within the same software ecosystem when transition from burst-neutrino calculations to DSNB predictions. My hope was that this would provide a common framework where users could: compute DSNB fluxes with the same physical conventions as SNEWPY, eventually use SNEWPY supernova emission models as DSNB source spectra directly rather than relying on fixed analytical spectra, compare burst-neutrino and diffuse-neutrino predictions within one package. If the maintainers think the current implementation is too independent from SNEWPY's existing architecture, I am happy to discuss redesigning it so that it interfaces directly with SNEWPY models, rather than existing as a largely standalone module.

[jpkneller]

Hi Dristii, Adding the capability to SNEWPY to compute the DSNB event rate in a detector is something that has been mentioned before so I am happy to see this contribution. I had a quick look through the code. My first impression is the same as Jost's: the overlap between dsnb.py and the rest of SNEWPY could be much greater than at the present time. There are two places in particular where I think the integration between them could be tighter: first in the supernova neutrino spectra, and second in using SNEWPY's RateCalculator (or it's interface with SNOwGLoBES) for the DSNB event rate calculation in more detectors and more detector channels. I have some thoughts on how to improve the integration but perhaps you have already better ideas.

[dristii-32]

Thank you for the feedback. I would very much like to hear your thoughts on how to improve the integration. The two directions you mentioned (tighter coupling with SNEWPY's emission models and with RateCalculator) are what I had in mind as next steps.

I am thinking that;

For the emission spectra: modify DSNBFlux to accept any SNEWPY supernova model (e.g. Nakazato_2013) as the sn_model argument, time-integrate its luminosity to obtain dN/dE per supernova, and use the current analytical PinchedSpectrum only as the default. This would mean a user could directly pass in a simulation-based spectrum rather than a parametric approximation.

For the event rate: replace the custom IBDCrossSection and integrated_ibd_rate with a call to SNEWPY's RateCalculator, inheriting support for all detector channels and materials already implemented there.

I can adapt the implementation around whichever interface fits best with SNEWPY's existing design. Any guidance on the intended architecture would be very helpful, and I’d be glad to work on a revised version based on your suggestions.

[jpkneller]

Those sound like the ideas I was going to suggest. You can find many scripts and notebooks in the doc directory that show how to use the supernova models and the RateCalculator. Here's one that might be useful https://github.com/SNEWS2/snewpy/blob/jpkneller-add_collate_for_RateCalculator/doc/source/nb/Analytic3Species.ipynb
A possible enhancement for version 2 of your code is to build a function that accepts a lot of supernova models and returns the average of their spectra according to some weight function e.g. the IMF. There are several model classes e.g. Warren_2020 where we have a lot of different simulation data sets: see https://github.com/SNEWS2/snewpy/blob/main/doc/source/nb/ccsn/Warren_2020.ipynb

[dristii-32]

Thank you for pointing to the Analytic3Species and Warren_2020 notebooks.

Based on your feedback, I am working on a revised version with two concrete improvements:

1. SNEWPY model integration for the source spectrum
   Rather than requiring users to specify a fixed analytical pinched spectrum, I will add a DSNBFlux.from_snewpy_model() classmethod that accepts any SNEWPY simulation model directly. Internally it time-integrates model.luminosity[flavor], model.meanE[flavor], and model.pinch[flavor] over the burst duration to obtain dN/dE per supernova, which is then used in the DSNB flux integral. This means a user can pass Nakazato_2013(...) or Warren_2020(...) directly as the source spectrum.

2. IMF-weighted model averaging
   I will also add DSNBFlux.from_snewpy_model_collection() which accepts a list of (model, progenitor_mass) pairs and computes the IMF-weighted average spectrum(the use case you described for the Warren_2020 model set where many progenitor masses are available). I plan to include a Salpeter IMF by default with an option to supply a custom weight function.

3. RateCalculator
   I would like to integrate with RateCalculator for the event rate step so users benefit from all detector channels already supported in SNEWPY.

Could you advise on the right interface, should DSNBFlux.integrated_ibd_rate() be replaced with a wrapper around RateCalculator, or would it be cleaner to keep the current analytic IBD rate as a fast default and offer RateCalculator as an optional backend?

I will push the revised code in the next few days.

[jpkneller]

1. and 2) sound good.
   WRT for 3), as shown in the Analytic3Species notebook, to use the RateCalculator you need to create a 'flux' or a 'fluence' object which you then pass as one of the arguments into RateCalculator's run function. A notebook which describes these objects can be found here https://github.com/SNEWS2/snewpy/blob/main/doc/source/nb/dev/FluxContainer_demo.ipynb and you might want to look into the source code which is in flux.py.

[JostMigenda]

1. SNEWPY model integration for the source spectrum
   Rather than requiring users to specify a fixed analytical pinched spectrum, I will add a DSNBFlux.from_snewpy_model() classmethod […] This means a user can pass Nakazato_2013(...) or Warren_2020(...) directly as the source spectrum.

What would be the use case for that?

We know that the DSNB spectrum will not be equivalent to the spectrum from any one progenitor; it is an IMF-weighted average over progenitor spectra with additional modifications (e.g., due to redshift). Existing spectra published in the literature take these into account to various degrees—they may still make unphysical assumptions, but at least those are described and motivated in the corresponding paper, and the limitations are generally well understood. Building an ad hoc spectrum based on an individual progenitor simulation lacks all that; so I don’t see a meaningful use case.

Without a clear use case, such a feature adds unnecessary complexity to the implementation and API—meaning more work for us and for users, as well as the potential for accidental misuse if inexperienced users don’t understand the difference between this and 2).

2. IMF-weighted model averaging
   I will also add DSNBFlux.from_snewpy_model_collection() which accepts a list of (model, progenitor_mass) pairs and computes the IMF-weighted average spectrum(the use case you described for the Warren_2020 model set where many progenitor masses are available). I plan to include a Salpeter IMF by default with an option to supply a custom weight function.

This sounds good to me!

Regarding the IMF: This is outside of my area of expertise, but I thought the differences between various IMFs are only in the low-mass region? If different IMFs all agree in the high-mass region relevant for SNe, do we actually need an option for custom IMFs? I’d probably prefer to drop this to keep the implementation and API simpler.

3. RateCalculator
   I would like to integrate with RateCalculator for the event rate step so users benefit from all detector channels already supported in SNEWPY.

Could you advise on the right interface, should DSNBFlux.integrated_ibd_rate() be replaced with a wrapper around RateCalculator, or would it be cleaner to keep the current analytic IBD rate as a fast default and offer RateCalculator as an optional backend?

I haven’t thought about the architecture yet; so I can’t comment on the design of the interface. But just as a general comment: I would like to avoid a situation where snewpy contains multiple different implementations of the same cross section. If the IBD implementation in RateCalculator has any significant issues, let’s tackle those directly, so that everyone benefits; rather than introduce a separate implementation that is only available for a fraction of overall functionality.

[dristii-32]

Thank you for the detailed feedback. I will revise the PR as follows:

Removing:

from_snewpy_model() for a single progenitor; agreed there is no physically meaningful use case for the DSNB.
IBDCrossSection class; agreed that a duplicate cross-section implementation is undesirable. Will remove entirely.
ibd_rate() and integrated_ibd_rate(); replaced by to_fluence() described below.

Simplifying:

from_snewpy_model_collection(); will retain this with a fixed Salpeter (1955) IMF and remove the custom IMF option to keep the API clean.

Adding:

DSNBFlux.to_fluence(energy, exposure); wraps the computed DSNB flux as a SNEWPY Fluence object that can be passed directly to RateCalculator.run(). This delegates all cross-section and detector calculations to the existing SNEWPY infrastructure, avoiding any duplicate implementations.

I will push the revised code shortly.

[dristii-32]

Hello, the revised code addresses all points raised:

IBDCrossSection and the duplicate cross-section implementation. ibd_rate() and integrated_ibd_rate() are gone. from_snewpy_model() for a single progenitor is removed.

from_snewpy_model_collection() now uses a fixed Salpeter (1955) IMF only. The custom IMF option is dropped.

DSNBFlux.to_fluence(energy, exposure) wraps the computed DSNB flux as a snewpy.flux.Fluence object. Users pass this directly to RateCalculator.run():

python
model = DSNBFlux()
fluence = model.to_fluence(np.linspace(10, 40, 200)u.MeV, exposure=10u.yr)
rates = RateCalculator().run(fluence)

This delegates all cross-section and detector calculations to the existing SNEWPY infrastructure with no duplicate implementations.

I wanted to ask that should DSNBFlux live at snewpy.dsnb (a new top-level module), or would you prefer it nested somewhere else in the package hierarchy, for example snewpy.models.dsnb? I will be glad to move it wherever fits best before this goes in.

[JostMigenda]

It’s late on Friday evening for me, so I haven’t looked in detail yet, but here are a couple of comments on my first scroll through. Many are fairly minor bits of cleanup, but there are a few open physics questions/issues that require further discussion.

[JostMigenda]

Suggested change

"IBDCrossSection",

Some cleanup. (Yep, forgetting to update __all__ happens to me regularly, too … 😅)

[JostMigenda]

I think these were only used for the IBD implementation?

[JostMigenda]

This looks like a leftover from the first, stand-alone version of the code. Can you clean this up?

(As well as at least one other instance of this, later in the file.)

[JostMigenda]

Also, since we know that snewpy.neutrino.Flavor exists, that export can move to the top of the file and we can set the default value of flavor in the definition right away, rather than checking for None in the function body.

[JostMigenda]

Suggested change

	>>> from snewpy_dsnb.dsnb import DSNBFlux
	>>> from snewpy.dsnb import DSNBFlux

[JostMigenda]

I think the weighting needs to be more complicated than this, in the general case. In this example, the 13 Msol progenitor would have a weight of 13^-2.35/(13^-2.35+20^-2.35) = 0.733. Now imagine that I accidentally add the 20 Msol progenitor to pairs twice—that would suddenly decrease the 13 Msol progenitor’s weight to 13^-2.35/(13^-2.35+2*20^-2.35) = 0.579.

(I’m wondering if a {mass: model} dictionary might be a slightly nicer data type here than a list of (model, mass) tuples? In part because it removes a level of nesting and in part because it eliminates that risk of accidental duplicates.)

Additionally, similar issues arise if I were to add a lot of progenitors in a small mass interval—that mass interval would suddenly have a much higher weight. Instead, I suspect the weight should be based on the integrated area under the IMF distribution that each progenitor represents? (Just my first guess; we should think this through more thoroughly than I can do on a Friday evening.)

[JostMigenda]

It looks like this is an accidental duplicate?

[JostMigenda]

Similar to the comment above: The exports can move to the top of the file and we can set the default value of flavor in the definition right away, rather than checking for None in the function body.

[JostMigenda]

This is no longer used, I think.

[JostMigenda]

This looks like a duplicate; class SNEWPYSpectrum is already defined further up in this file? Though I’m not sure if there are differences between the two (and if so, which is the correct one).

[dristii-32]

Hello, all 21 tests now pass. Summary of changes:

• from_snewpy_model (single progenitor) removed

• IBDCrossSection and all duplicate cross-section code removed

• ibd_rate() and integrated_ibd_rate() removed; replaced by to_fluence()

• from_snewpy_model_collection now accepts a {mass: model} dict; uses Salpeter IMF integrated over proper mass intervals (midpoints between adjacent progenitor masses) rather than point evaluation: physically correct for uneven mass spacing
  On f_bh in from_snewpy_model_collection: you are right. When the collection already contains both successful and BH-forming progenitor models, the BH fraction is implicit in the set of models provided and their IMF weights — a separate f_bh parameter is redundant and potentially misleading. Removed it from the method signature; from_snewpy_model_collection now constructs the full IMF-weighted spectrum directly from whatever models are passed in.

• All try/except ImportError blocks for snewpy removed; Flavor moved to top-level import

• All example imports updated to from snewpy.dsnb import DSNBFlux

• Test suite updated to match the revised API

[dristii-32]

Hi, just checking in.. Happy to make any further changes if needed... Let me know if there is anything else before this is ready for merge.

[jpkneller]

@dristii-32, Sorry for the delay in getting back to you.

My personal preference is that the DSNB flux becomes a model class similar to, say, Bollig_2016 or Warren_2020. Method names for the DSNB model should be identical to those for all the other models (in those cases where it makes sense to have similar methods) so that the model is intuitive for users. We definitely need a usage example which might be along the lines of:

    model = DSNB(Warren_2020) # use the Warren_2020 models to build the DSNB 
    transformation = AdiabaticMSW(MixingParameters('NORMAL')) # Desired flavor transformation
   
    times    = np.linspace(0,1,2) << u.yr # beginning and end of detector exposure
    energies = np.linspace(0,50,501) << u.MeV

    #get the flux from the model
    flux = model.get_flux(t=times, E=energies, flavor_xform=transformation)
    fluence = flux.integrate('time') 

    detector = "wc100kt30prct"
    rc = RateCalculator()

    events = rc.run(fluence, detector)        

Perhaps @JostMigenda has ideas too. The model = DSNB(Warren_2020) is meant to be suggestive: more parameters might need to be supplied or maybe the input to the class is a list of supernova model instances the user has previously created along with a way to weight them. Is that clear?

Perhaps I don't fully understand the code but as Jost said before, when I look through dsnb.py I can see a lot of code which looks to a) duplicate existing parts of SNEWPY and b) does not integrate well with SNEWPY because it does not use the essential components such as the Container for spectra / fluxes / fluences. In particular, I am confused by the need for the SNEWPYSpectrum and PinchedSpectrum classes. SNEWPYSpectrum appears to time integrate the spectrum (shouldn't it time average?) but only for a single flavor and without a flavor transformation prescription. It also only works for pinched models. This same functionality (but without the limitations) can be achieved with spectra = model.get_flux and spectra.integrate('time') as in the example above. PinchedSpectrum is also a very similar to SNEWPY's existing model class PinchedModel.

[dristii-32]

Thank you for the redesign direction. I understand that I should

• Rename the class to DSNB and make get_flux(t, E, flavor_xform) the primary method, matching the interface of other SNEWPY model classes
• Remove PinchedSpectrum, SNEWPYSpectrum, and CoreCollapseRate as standalone classes; the SN emission spectrum comes directly from SNEWPY model instances via their existing interface
• The t argument to get_flux is the detector exposure window; the DSNB flux is steady-state so flux.integrate('time') over [0, T] gives the fluence over exposure time T
• And apply flavor_xform to the SN emission spectrum before the redshift integration

Before I implement, I wanted to ask that for DSNB(Warren_2020) where Warren_2020 is the class; should the code automatically enumerate all available progenitor parameter combinations from Warren_2020.param and apply IMF weighting by progenitor mass...? Or is the preferred interface DSNB([instance1, instance2, ...]) where the user creates the instances themselves?

The first is more convenient but requires knowing which parameter key corresponds to progenitor mass (which varies across model classes). The second is more explicit. Or maybe I can also implement a combination where passing a class triggers the auto-enumeration and passing a list uses those instances directly.

I will begin the rewrite once this is clear.

[jpkneller]

* Remove PinchedSpectrum, SNEWPYSpectrum, and CoreCollapseRate as standalone classes; the SN emission spectrum comes directly from SNEWPY model instances via their existing interface

From what I can see, PinchedSpectrum and SNEWPYSpectrum duplicate existing code or existing methods. CoreCollapseRate can be kept.

* The t argument to get_flux is the detector exposure window; the DSNB flux is steady-state so flux.integrate('time') over [0, T] gives the fluence over exposure time T

Yes, I would suggest that option so that the interface with the DSNB model is most similar to any supernova model.

Before I implement, I wanted to ask that for DSNB(Warren_2020) where Warren_2020 is the class; should the code automatically enumerate all available progenitor parameter combinations from Warren_2020.param and apply IMF weighting by progenitor mass...? Or is the preferred interface DSNB([instance1, instance2, ...]) where the user creates the instances themselves?
The first is more convenient but requires knowing which parameter key corresponds to progenitor mass (which varies across model classes). The second is more explicit. Or maybe I can also implement a combination where passing a class triggers the auto-enumeration and passing a list uses those instances directly.

I agree with you: my preference is also for the second option because it gives the user the greatest flexibility but I can see that it might make the code more complicated. There is still the potential for a problem if more than one model in the list has the same mass.

I was not aware that the progenitor mass key was different across the models: we will fix that. For now use the key for the Warren_2020 models.

[dristii-32]

Hello, the revised code has been pushed. The module now has a single DSNB class with get_flux(t, E, flavor_xform) as the primary method, matching the interface of other SNEWPY model classes. PinchedSpectrum and SNEWPYSpectrum are removed; the per-supernova emission spectrum is computed directly from each SNEWPY model's luminosity, meanE, and pinch time series via time integration. CoreCollapseRate is kept.

Usage now follows this pattern:

python
model = DSNB([(sn1, 13), (sn2, 20)])
flux = model.get_flux(t=[0,1]*u.yr, E=energies,
flavor_xform=NoTransformation())
fluence = flux.integrate('time')
events = RateCalculator().run(fluence, "wc100kt30prct")

Also, flavor_xform is currently reserved for future use. The NU_E_BAR component carries the full DSNB flux since MSW effects average out over cosmological distances for the standard calculation.

[jpkneller]

Great! I'll take a look.
I think there may be some confusion over the flavor transformation that you need to provide in get_flux. What this flavor transformation describes is the effect of the supernova matter upon the neutrinos. This does not average out - every DSNB calculation e.g. https://arxiv.org/pdf/1205.6292 distinguishes between the two mass orderings.

[dristii-32]

Thank you for the pointer. The revised code now,

• Applies flavor_xform inside _dNdE_per_sn, before the redshift integration, using the oscillated spectrum dN/dE_osc = p_ee · dN/dE_nuebar + p_xe · dN/dE_nux
• Uses the correct SNEWPY API: prob_eebar(t, E) and prob_xebar(t, E) from the flavor transformation object
• Pass AdiabaticMSW(mh=MassHierarchy.NORMAL) or AdiabaticMSW(mh=MassHierarchy.INVERTED) for the two hierarchies; NoTransformation() or None gives the unoscillated flux

Also, AdiabaticMSW.prob_eebar returns a scalar independent of (t, E), so _oscillation_probs queries it once at a reference point. If a time/energy-dependent transformation is needed in the future, _dNdE_per_sn would need to evaluate probabilities on the full (N_t, N_E) grid. I will be happy to add that if required.

Tests confirm NH/IH ratio and the ~50–60% MSW effect on the NU_E_BAR flux, consistent with Lunardini & Tamborra (2012).