Fornax 2021 Models

Neutrino spectra from the long-time run 2D (axisymmetric) models produced by Burrows and Vartanyan, Nature 589:29-39, 2021.

Data taken from the HDF5 files available for download at the Princeton group website.

from snewpy.neutrino import Flavor
from snewpy.models.ccsn import Fornax_2021

from astropy import units as u
from glob import glob

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

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

Initialize Models

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

Fornax_2021.param

{'progenitor_mass': <Quantity [12.  , 13.  , 14.  , 15.  , 16.  , 17.  , 18.  , 19.  , 20.  ,
            21.  , 22.  , 23.  , 25.  , 26.  , 26.99] solMass>}

We’ll initialise some of these progenitors and plot the luminosity of different neutrino flavors for two of them. (Note that the Fornax_2021 simulations didn’t distinguish between \(\nu_x\) and \(\bar{\nu}_x\), so both have the same luminosity.) If this is the first time you’re using a progenitor, snewpy will automatically download the required data files for you.

models = {}
for m in Fornax_2021.param['progenitor_mass'].value[::2]:
    # Initialise every second progenitor
    models[m] = Fornax_2021(progenitor_mass=m*u.solMass)

models[12]

Fornax_2021 Model

Parameter	Value
Progenitor mass	\(12\) \(\mathrm{M_{\odot}}\)

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, 1.0),
           xlabel=r'$t-t_{\rm bounce}$ [s]',
           title=r'{} $M_\odot$'.format(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}$]');

Spectra of All Flavors vs. Time for the \(12M_{\odot}\) Model

Use Default Linear Interpolation in Flux Retrieval

model = models[12]  # Use the 12 solar mass model

times = np.arange(-0.2, 3.8, 0.2) * u.s
E = np.arange(0, 101, 1) * u.MeV

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

linestyles = ['-', '--', '-.', ':']

spectra = model.get_initial_spectra(times, E)

for i, ax in enumerate(axes.flatten()):
    for line, flavor in zip(linestyles, Flavor):
        ax.plot(E, spectra[flavor][i], lw=3, ls=line, label=flavor.to_tex())
    ax.set(xlim=(0,100))
    ax.set_title('$t$ = {:g}'.format(times[i]), fontsize=16)
    ax.legend(loc='upper right', ncol=2, fontsize=12)

fig.text(0.5, 0., 'energy [MeV]', ha='center')
fig.text(0., 0.5, f'flux [{spectra[Flavor.NU_E].unit}]', va='center', rotation='vertical');

Use Nearest-Bin “Interpolation” in Flux Retrieval

times = np.arange(-0.2, 3.8, 0.2) * u.s
E = np.arange(0, 101, 1) * u.MeV

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

linestyles = ['-', '--', '-.', ':']

spectra = model.get_initial_spectra(times, E, interpolation='nearest')

for i, ax in enumerate(axes.flatten()):
    for line, flavor in zip(linestyles, Flavor):
        ax.plot(E, spectra[flavor][i], lw=3, ls=line, label=flavor.to_tex())
    ax.set(xlim=(0,100))
    ax.set_title('$t$ = {:g}'.format(times[i]), fontsize=16)
    ax.legend(loc='upper right', ncol=2, fontsize=12)

fig.text(0.5, 0., 'energy [MeV]', ha='center')
fig.text(0., 0.5, f'flux [{spectra[Flavor.NU_E].unit}]', va='center', rotation='vertical');

Progenitor Mass Dependence

Luminosity vs. Time for a Selected List of Progenitor Masses

Plot \(L_{\nu_{e}}(t)\) for a selection of progenitor masses to observe the dependence of the emission on mass.

fig, axes = plt.subplots(3,1, figsize=(10,13), sharex=True, sharey=True,
                         gridspec_kw = {'hspace':0.02})

colors0 = mpl.cm.viridis(np.linspace(0.1,0.9, len(models)))
colors1 = mpl.cm.inferno(np.linspace(0.1,0.9, len(models)))
colors2 = mpl.cm.cividis(np.linspace(0.1,0.9, len(models)))

linestyles = ['-', '--', '-.', ':']

for i, model in enumerate(models.values()):
    ax = axes[0]
    flavor = Flavor.NU_E
    ax.plot(model.time, model.luminosity[flavor], lw=2, color=colors0[i], ls=linestyles[i%4],
            label='${0.value:g}$ {0.unit:latex}'.format(model.progenitor_mass))
    ax.set(xscale='log',
           xlim=(1e-3, 4),
           yscale='log',
           ylim=(0.4e52, 9e53),
           ylabel=r'$L_{\nu_e}(t)$ [erg s$^{-1}$]')
    ax.grid(ls=':', which='both')
    ax.legend(ncol=3, fontsize=12, title=r'$\nu_e$');

    ax = axes[1]
    flavor = Flavor.NU_E_BAR
    ax.plot(model.time, model.luminosity[flavor], lw=2, color=colors1[i], ls=linestyles[i%4],
        label='${0.value:g}$ {0.unit:latex}'.format(model.progenitor_mass))
    ax.set(ylabel=r'$L_{\bar{\nu}_e}(t)$ [erg s$^{-1}$]')
    ax.grid(ls=':', which='both')
    ax.legend(ncol=3, fontsize=12, title=r'$\bar{\nu}_e$');

    ax = axes[2]
    flavor = Flavor.NU_X
    ax.plot(model.time, model.luminosity[flavor], lw=2, color=colors2[i], ls=linestyles[i%4],
        label='${0.value:g}$ {0.unit:latex}'.format(model.progenitor_mass))
    ax.set(xlabel='time [s]',
           ylabel=r'$L_{\nu_X}(t)$ [erg s$^{-1}$]')
    ax.grid(ls=':', which='both')
    ax.legend(ncol=3, fontsize=12, title=r'$\nu_X$');

Progenitor Dependence of Spectra at 70 ms

Use Default Linear Interpolation in Flux Retrieval

t = 70*u.ms
E = np.arange(0, 101, 1) * u.MeV

fig, axes = plt.subplots(2,4, figsize=(16,6), sharex=True, sharey=True, tight_layout=True)

linestyles = ['-', '--', '-.', ':']

for model, ax in zip(models.values(), axes.flatten()):
    spectra = model.get_initial_spectra(t, E)
    for line, flavor in zip(linestyles, Flavor):
        ax.plot(E, spectra[flavor][0], lw=3, ls=line, label=flavor.to_tex())
    ax.set(xlim=(0,100))
    ax.set_title('${0.value:g}$ {0.unit:latex}'.format(model.progenitor_mass))
    ax.legend(loc='upper right', ncol=2, fontsize=12)
    ax.grid(ls=':')

fig.text(0.5, 0., 'energy [MeV]', ha='center')
fig.text(0., 0.5, f'flux [{spectra[Flavor.NU_E].unit}]', va='center', rotation='vertical');

Use Nearest-Bin “Interpolation” in Flux Retrieval

t = 70*u.ms
E = np.arange(0, 101, 1) * u.MeV

fig, axes = plt.subplots(2,4, figsize=(16,6), sharex=True, sharey=True, tight_layout=True)

linestyles = ['-', '--', '-.', ':']

for model, ax in zip(models.values(), axes.flatten()):
    spectra = model.get_initial_spectra(t, E, interpolation='nearest')
    for line, flavor in zip(linestyles, Flavor):
        ax.plot(E, spectra[flavor][0], lw=3, ls=line, label=flavor.to_tex())
    ax.set(xlim=(0,100))
    ax.set_title('${0.value:g}$ {0.unit:latex}'.format(model.progenitor_mass))
    ax.legend(loc='upper right', ncol=2, fontsize=12)
    ax.grid(ls=':')

fig.text(0.5, 0., 'energy [MeV]', ha='center')
fig.text(0., 0.5, f'flux [{spectra[Flavor.NU_E].unit}]', va='center', rotation='vertical');