Flavor Transformations¶
This notebook produces the figures in the SNEWPY paper showing the effect of the flavor transformation prescriptions upon the nakazato-shen-z0.004-t_rev100ms-s20.0 model
[1]:
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, MixingParameters
from snewpy.models.ccsn import Nakazato_2013
from snewpy.flavor_transformation import AdiabaticMSW, NonAdiabaticMSWH, \
TwoFlavorDecoherence, ThreeFlavorDecoherence, \
NeutrinoDecay, AdiabaticMSWes, NonAdiabaticMSWes
mpl.rc('font', size=18)
%matplotlib inline
# `np.trapz` was renamed in numpy 2.0 and now raises a DeprecationWarning
if np.__version__ < "2.0":
np.trapezoid = np.trapz
[2]:
model = Nakazato_2013(progenitor_mass=20*u.solMass, revival_time=100*u.ms, metallicity=0.004, eos='shen')
model
[2]:
Nakazato_2013 Model: nakazato-shen-z0.004-t_rev100ms-s20.0.fits
Parameter |
Value |
|---|---|
Progenitor mass |
\(20\) \(\mathrm{M_{\odot}}\) |
Revival time |
\(100\) \(\mathrm{ms}\) |
Metallicity |
0.004 |
EOS |
shen |
Transformed and untransformed energy integrated flux at Earth¶
Compute and plot the flux at Earth for a SN at 10 kpc with no flavor transformation, and with the chosen flavor transformation.
[3]:
def plot_total_flux(model, xform_nmo, xform_imo):
"""Plot initial and oscillated neutrino luminosities.
Parameters
----------
model : SupernovaModel
An input model from a CCSN simulation.
flav_xform : FlavorTransformation
A FlavorTransformation subclass; used to create an instance.
"""
energies = np.linspace(0,60,121) * u.MeV
d = (10*u.kpc).to('cm').value # distance to SN
times = model.get_time()
burst_epoch = times <= 0.1*u.s
accretion_epoch = (times > 0.1*u.s) & (times <= 0.5*u.s)
cooling_epoch = (times > 0.5*u.s) & (times <= 10*u.s)
ilum = {}
olum_nmo = {}
olum_imo = {}
for flavor in Flavor:
ilum[flavor] = np.zeros(len(times))
olum_nmo[flavor] = np.zeros(len(times))
olum_imo[flavor] = np.zeros(len(times))
# Compute the transformed and untransformed flux at each time.
for i, t in enumerate(times):
ispec = model.get_initial_spectra(t, energies)
ospec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)
ospec_imo = model.get_transformed_spectra(t, energies, xform_imo)
for flavor in Flavor:
for j, E in enumerate(energies):
ispec[flavor][j] /= (4.*np.pi*d**2)
ospec_nmo[flavor][j] /= (4.*np.pi*d**2)
ospec_imo[flavor][j] /= (4.*np.pi*d**2)
for flavor in Flavor:
ilum[flavor][i] = np.trapezoid(ispec[flavor].to('1/(erg*s)'), energies.to('erg')).value
olum_nmo[flavor][i] = np.trapezoid(ospec_nmo[flavor].to('1/(erg*s)'), energies.to('erg')).value
olum_imo[flavor][i] = np.trapezoid(ospec_imo[flavor].to('1/(erg*s)'), energies.to('erg')).value
# make the figures
fig, axes = plt.subplots(3,3, figsize=(20,12), tight_layout=True)
smax = [0.,0.,0.]
titles = ['Untransformed', 'Transformed (NMO)', 'Transformed (IMO)']
for i, spec in enumerate([ilum, olum_nmo, olum_imo]):
for j, phase in enumerate([burst_epoch, accretion_epoch, cooling_epoch]):
ax = axes[i,j]
timeunits = 'ms' if j==0 else 's'
for flavor in Flavor:
if i == 0:
smax[j] = np.maximum(smax[j], 1.1*np.max(spec[flavor][phase]))
ax.plot(times[phase].to(timeunits),
spec[flavor][phase], label=flavor.to_tex(), lw=3,
color='C0' if flavor.is_electron else 'C1',
ls='-' if flavor.is_neutrino else ':')
ax.set(xlim=(times[phase][0].to(timeunits).value, times[phase][-1].to(timeunits).value),
ylim=(0, smax[j]))
if j==0:
ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')
ax.legend(loc='upper right', ncol=1, fontsize=18)
if j==1:
ax.set(title=titles[i])
if i < 2:
ax.set(xticklabels=[])
else:
ax.set(xlabel='time [{}]'.format(timeunits))
ax.grid(ls=':')
return fig
Untransformed and Transformed Spectra at Earth¶
Compute and plot the spectra at Earth for a SN at 10 kpc with no flavor transformation, and with the chosen flavor transformation.
[4]:
def plot_spectra(model, xform_nmo, xform_imo, t):
"""Plot initial and oscillated neutrino luminosities.
Parameters
----------
model : SupernovaModel
An input model from a CCSN simulation.
flav_xform : FlavorTransformation
A FlavorTransformation subclass; used to create an instance.
t : astropy.Quantity
Time to compute the spectrum.
"""
energies = np.linspace(0,60,121) * u.MeV
d = (10*u.kpc).to('cm').value # distance to SN
#get the spectra
ispec = model.get_initial_spectra(t, energies)
ospec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)
ospec_imo = model.get_transformed_spectra(t, energies, xform_imo)
for flavor in Flavor:
for j, E in enumerate(energies):
ispec[flavor][j] /= (4.*np.pi*d**2)
ospec_nmo[flavor][j] /= (4.*np.pi*d**2)
ospec_imo[flavor][j] /= (4.*np.pi*d**2)
fig, axes = plt.subplots(2,2, figsize=(18,15), sharex=True, sharey=True, tight_layout=True)
for i, spec in enumerate([ispec, ospec_nmo]):
axes[0][0].plot(energies, spec[Flavor.NU_E]/1e16,
label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),
color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)
axes[0][0].plot(energies, spec[Flavor.NU_X]/1e16,
label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),
color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)
axes[0][0].set(title='Neutrinos in the NMO: $t = ${:.1f}'.format(t))
axes[0][0].grid()
axes[0][0].legend(loc='upper right', ncol=2, fontsize=16)
axes[0][1].plot(energies, spec[Flavor.NU_E_BAR]/1e16,
label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),
color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)
axes[0][1].plot(energies, spec[Flavor.NU_X_BAR]/1e16,
label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),
color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)
axes[0][1].set(title='Antineutrinos in the NMO: $t = ${:.1f}'.format(t))
axes[0][1].grid()
axes[0][1].legend(loc='upper right', ncol=2, fontsize=16)
for i, spec in enumerate([ispec, ospec_imo]):
axes[1][0].plot(energies, spec[Flavor.NU_E]/1e16,
label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),
color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)
axes[1][0].plot(energies, spec[Flavor.NU_X]/1e16,
label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),
color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)
axes[1][0].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Neutrinos in the IMO: $t = ${:.1f}'.format(t))
axes[1][0].grid()
axes[1][0].legend(loc='upper right', ncol=2, fontsize=16)
axes[1][1].plot(energies, spec[Flavor.NU_E_BAR]/1e16,
label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),
color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)
axes[1][1].plot(energies, spec[Flavor.NU_X_BAR]/1e16,
label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),
color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)
axes[1][1].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Antineutrinos in the IMO: $t = ${:.1f}'.format(t))
axes[1][1].grid()
axes[1][1].legend(loc='upper right', ncol=2, fontsize=18)
ax = axes[0][0]
ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')
ax = axes[1][0]
ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')
return fig
[5]:
fig = plot_total_flux(model, AdiabaticMSW(), AdiabaticMSW(mh=MassHierarchy.INVERTED) )
fig.show()
# fig.savefig('flux_adiabaticmsw.pdf')
fig = plot_spectra(model, AdiabaticMSW(), AdiabaticMSW(mh=MassHierarchy.INVERTED), 100*u.ms)
fig.show()
# fig.savefig('spectra_adiabaticmsw.pdf')
[6]:
mp_nmo = MixingParameters()
angles_nmo = mp_nmo.get_mixing_angles() + (10*u.deg,)
xf_nmo = AdiabaticMSWes(angles_nmo)
mp_imo = MixingParameters(mh=MassHierarchy.INVERTED)
angles_imo = mp_imo.get_mixing_angles() + (10*u.deg,)
xf_imo = AdiabaticMSWes(angles_imo)
fig = plot_total_flux(model, xf_nmo, xf_imo)
fig.show()
# fig.savefig('flux_adiabaticmswes.pdf')
fig = plot_spectra(model, xf_nmo, xf_imo, 100*u.ms)
fig.show()
# fig.savefig('spectra_adiabaticmswes.pdf')
