O’Connor 2013 Models

Data from O’Connor & Ott 2013, 32 progenitors (Woosley and Heger 2007) and 2 EOS (LS220 and HShen) for 500 ms post bounce in spherical symmetry (no explosions)

Reference: O’Connor and Ott ApJ 762 126 2013

• doi:10.1088/0004-637X/762/2/126

• arXiv:1207.1100

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 OConnor_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 OConnor_2013 model. We can use the param property to view all physics parameters and their possible values:

OConnor_2013.param

{'progenitor_mass': <Quantity [ 12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,
             22.,  23.,  24.,  25.,  26.,  27.,  28.,  29.,  30.,  31.,
             32.,  33.,  35.,  40.,  45.,  50.,  55.,  60.,  70.,  80.,
            100., 120.] solMass>,
 'eos': ['HShen', 'LS220']}

Quite a lot of choice there! Let’s initialise all of these progenitors and compare their \(\nu_e\) luminosities. If this is the first time you’re using a progenitor, snewpy will automatically download the required data files for you.

models = {}
for mass in OConnor_2013.param['progenitor_mass']:
    models[int(mass.value)] = OConnor_2013(progenitor_mass=mass, eos='LS220')

for model in models.values():
    plt.plot(model.time, model.luminosity[Flavor.NU_E]/1e51, 'C0', lw=1)

plt.xlabel(r'$t$ [s]')
plt.ylabel(r'luminosity [foe s$^{-1}$]');

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

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

for i, model in enumerate([models[12], models[20]]):
    ax = axes[i]
    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.51),
           xlabel=r'$t-t_{\rm bounce}$ [s]',
           title=r'{}: {} $M_\odot$'.format(model.metadata['EOS'], model.metadata['Progenitor mass'].value))
    ax.grid()
    ax.legend(loc='upper right', ncol=2, fontsize=18)

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

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 = 400*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}$]')

fig.tight_layout();