Nakazato 2013 Models

This notebook includes an example of how to read the .fits files located in same directory as this notebook. These files contain neutrino luminosity, mean energy, and energy spectrum pinch parameter (alpha) as functions of time for the CCSN neutrino flavors. These were obtained using models from Nakazato et al., 2013 and 2015. Data are public and taken from their website.

The citation for use of the database is: Supernova Neutrino Light Curves and Spectra for Various Progenitor Stars: From Core Collapse to Proto-neutron Star Cooling, K. Nakazato, K. Sumiyoshi, H. Suzuki, T. Totani, H. Umeda, and S. Yamada, Astrophys. J. Supp. 205 (2013) 2, arXiv:1210.6841.

If the BH model with LS220 EOS is used, the citation is: Spectrum of the Supernova Relic Neutrino Background and Metallicity Evolution of Galaxies, K. Nakazato, E. Mochida, Y. Niino, and H. Suzuki, Astrophys. J. 804 (2015) 75, arXiv:1503.01236.

For the BH model with Togashi EOS, the citation is: Numerical Study of Stellar Core Collapse and Neutrino Emission Using the Nuclear Equation of State Obtained by the Variational Method, K. Nakazato, K. Sumiyoshi, and H. Togashi, Publ. Astron. Soc. Jpn. 73 (2021) 639-651, arXiv:2103.14386.

These examples use code taken from IceCube’s fast supernova monte carlo, ASTERIA.

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

from astropy import units as u

from snewpy.neutrino import Flavor, MassHierarchy
from snewpy.models.ccsn import Nakazato_2013
from snewpy.flavor_transformation import NoTransformation, AdiabaticMSW, ThreeFlavorDecoherence

mpl.rc('font', size=16)
%matplotlib inline

Initialize Models

To start, let’s see what progenitors are available for the Nakazato_2013 model. We can use the param property to view all physics parameters and their possible values:

Nakazato_2013.param

{'progenitor_mass': <Quantity [13., 20., 30., 50.] solMass>,
 'revival_time': <Quantity [  0., 100., 200., 300.] ms>,
 'metallicity': [0.02, 0.004],
 'eos': ['LS220', 'shen', 'togashi']}

Quite a lot of choice there! However, for the Nakazato model, not all combinations of these parameters are valid. We can use the get_param_combinations function to get a list of all valid combinations or to filter it:

# This will print a long tuple of combinations:
#Nakazato_2013.get_param_combinations()

# To get a more manageable list, let’s filter it to show only the
# available progenitors with a metallicity of z=0.004:
highz_models = list(params for params in Nakazato_2013.get_param_combinations() if params['metallicity'] == 0.02)
print("Progenitors with metallicity z=0.02:", *highz_models, sep='\n')

Progenitors with metallicity z=0.02:
{'progenitor_mass': <Quantity 13. solMass>, 'revival_time': <Quantity 100. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 13. solMass>, 'revival_time': <Quantity 200. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 13. solMass>, 'revival_time': <Quantity 300. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 20. solMass>, 'revival_time': <Quantity 100. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 20. solMass>, 'revival_time': <Quantity 200. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 20. solMass>, 'revival_time': <Quantity 300. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 30. solMass>, 'revival_time': <Quantity 100. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 30. solMass>, 'revival_time': <Quantity 200. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 30. solMass>, 'revival_time': <Quantity 300. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 50. solMass>, 'revival_time': <Quantity 100. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 50. solMass>, 'revival_time': <Quantity 200. ms>, 'metallicity': 0.02, 'eos': 'shen'}
{'progenitor_mass': <Quantity 50. solMass>, 'revival_time': <Quantity 300. ms>, 'metallicity': 0.02, 'eos': 'shen'}

We’ll pick one of these progenitors and initialise it. If this is the first time you’re using this progenitor, snewpy will automatically download the required data file for you.

model_params = highz_models[3]
model = Nakazato_2013(**model_params)
# This is equivalent to:
#model = Nakazato_2013(progenitor_mass=20*u.solMass, revival_time=100*u.ms, metallicity=0.004, eos='shen')
model

Nakazato_2013 Model: nakazato-shen-z0.02-t_rev100ms-s20.0.fits

Parameter	Value
Progenitor mass	\(20\) \(\mathrm{M_{\odot}}\)
Revival time	\(100\) \(\mathrm{ms}\)
Metallicity	0.02
EOS	shen

Finally, let’s plot the luminosity of different neutrino flavors for this model. (Note that the Nakazato_2013 simulations didn’t distinguish between \(\nu_x\) and \(\bar{\nu}_x\), so both have the same luminosity.)

fig, ax = plt.subplots(1, figsize=(8,6), tight_layout=False)

for flavor in Flavor:
    ax.plot(model.time, model.luminosity[flavor]/1e51, # Report luminosity in units foe/s
            label=flavor.to_tex(),
            color = 'C0' if flavor.is_electron else 'C1',
            ls = '-' if flavor.is_neutrino else ':',
            lw = 2 )

ax.set(xlim=(-0.05, 0.5),
       xlabel=r'$t-t_{\rm bounce}$ [s]',
       ylabel=r'luminosity [foe s$^{-1}$]')
ax.grid()
ax.legend(loc='upper right', ncol=2, fontsize=18);

Initial and Oscillated Spectra

Plot the neutrino spectra at the source and after the requested flavor transformation has been applied.

Adiabatic MSW Flavor Transformation: Normal mass ordering

# Adiabatic MSW effect. NMO is used by default.
xform_nmo = AdiabaticMSW()

# Energy array and time to compute spectra.
# Note that any convenient units can be used and the calculation will remain internally consistent.
E = np.linspace(0,100,201) * u.MeV
t = 100*u.ms

ispec = model.get_initial_spectra(t, E)
ospec_nmo = model.get_transformed_spectra(t, E, xform_nmo)

fig, axes = plt.subplots(1,2, figsize=(12,5), sharex=True, sharey=True, tight_layout=True)

for i, spec in enumerate([ispec, ospec_nmo]):
    ax = axes[i]
    for flavor in Flavor:
        ax.plot(E, spec[flavor],
                label=flavor.to_tex(),
                color='C0' if flavor.is_electron else 'C1',
                ls='-' if flavor.is_neutrino else ':', lw=2,
                alpha=0.7)

    ax.set(xlabel=r'$E$ [{}]'.format(E.unit),
           title='Initial Spectra: $t = ${:.1f}'.format(t) if i==0 else 'Oscillated Spectra: $t = ${:.1f}'.format(t))
    ax.grid()
    ax.legend(loc='upper right', ncol=2, fontsize=16)

ax = axes[0]
ax.set(ylabel=r'flux [erg$^{-1}$ s$^{-1}$]');