Fornax 2019 Models

Some I/O and plotting using the Fornax 3D models from Vartanyan, Burrows, et al., MNRAS 482(1):351, 2019, which express the angular dependence of the neutrino luminosity from CCSNe in terms of a spherical harmonic expansion up to degree \(\ell=2\). Note: Caching and plotting the angular dependence requires the healpy package, which can be installed via pip install healpy.

The data are available on the Burrows group website.

[1]:

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

from astropy import units as u

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

try:
    import healpy as hp
    healpy_available = True
except ModuleNotFoundError:
    healpy_available = False

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

Initialization and Time Evolution of the Energy Spectrum

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

[2]:

Fornax_2019.param

[2]:

{'progenitor_mass': <Quantity [ 9., 10., 12., 13., 14., 15., 16., 19., 25., 60.] solMass>}

We’ll initialise one of these progenitors and plot the evolution of the predicted neutrino spectrum across four time steps. If this is the first time you’re using a progenitor, snewpy will automatically download the required data files for you.

[3]:

model = Fornax_2019(progenitor_mass=10*u.solMass)
model

[3]:

Fornax_2019 Model: lum_spec_10M.h5

Parameter

Value

Progenitor mass

\(10\) \(\mathrm{M_{\odot}}\)

 
[4]:

E = np.arange(0,111) * u.MeV
theta = 23.55*u.degree
phi = 22.5*u.degree

fig, axes = plt.subplots(2,2, figsize=(11,8), sharex=True, sharey=True,
                         gridspec_kw = {'hspace':0.04, 'wspace':0.025})
axes = axes.flatten()

for ax, t in zip(axes, [10., 100., 300., 600.]):
    t *= u.ms
    spectra = model.get_initial_spectra(t, E, theta, phi)

    lines = ['--', '-.', '-', ':']
    for line, flavor in zip(lines, Flavor):
        ax.plot(E, spectra[flavor].to('1e52 erg/(s MeV)'), ls=line, label=flavor.to_tex())
    ax.legend(fontsize=14, ncol=2, loc='upper right', title='{:g}'.format(t))
    ax.set(xlim=(E[0].to_value('MeV'), E[-1].to_value('MeV')))
    ax.grid(ls=':')

axes[0].set(ylabel=r'$dL/dE$ [$10^{52}$ erg s$^{-1}$ MeV$^{-1}$]')
axes[2].set(xlabel=r'energy [MeV]');

fig.suptitle('Time Evolution of Spectra'.format(t))
fig.subplots_adjust(top=0.925, bottom=0.1)
../../_images/nb_ccsn_Fornax_2019_6_0.png

Test Caching of Full Model Angular Dependence

Here, if the user specifies cache_flux=True in the model constructor, the entire angular dependence of the model will be computed using a HEALPix grid with 192 pixels, using a HEALPix grid with parameter nside=4 and ordering=RING. The model constructor will calculate the CCSN fluxes on the unit sphere at the centers of the 192 healpixels for all energy and time bins in the model.

The resulting tables are then saved to a FITS file in the same location as the input HDF5 file. If the user instantiates this model with cache_flux=True in the future, the model will be initialized using the cached FITS file.

Comments:

  1. The HDF5 files store the angular dependence of the models in a highly compressed format using the spherical harmonic expansion up to degree \(\ell=2\).

  2. The time bins are equal across all flavors but the energy binning is a function of both time and flavor.

  3. Computing the fluxes at 192 locations on the sky takes about 30 seconds per flavor due to the overhead of reading out the HDF5 files. So construction of the cached model is slow when called for the first time.

  4. With nside=4 and 192 bins, the angular resolution per bin is about \(15^\circ\), which is reasonable given the highest-resolution features of the models are the quadrupole moments (\(\ell=2\)).

[5]:

cmodel = Fornax_2019(progenitor_mass=10*u.solMass, cache_flux=True)
cmodel

[5]:

Fornax_2019 Model: lum_spec_10M.h5

Parameter

Value

Progenitor mass

\(10\) \(\mathrm{M_{\odot}}\)

Test Equivalence of the Cached and Uncached Models

[6]:

t = 100 * u.ms

Es, dEs, spectras = model._get_binnedspectra(t, theta, phi)
Ec, dEc, spectrac = cmodel._get_binnedspectra(t, theta, phi)

for flavor in Flavor:
    print('\n{}'.format(str(flavor)))
    print(Es[flavor])
    print(Ec[flavor])
    print(spectras[flavor])
    print(spectrac[flavor])

    j = np.argmax(spectrac[flavor])
    print('Max value:')
    print(Ec[flavor][j], spectrac[flavor][j])


0
[  1.12409548   1.80813134   2.90841747   4.67825097   7.52506555
  12.10422696  19.46990488  31.31775346  50.37526829  81.02968365
 130.33795859 209.65136088] MeV
[  1.12409548   1.80813134   2.90841747   4.67825097   7.52506555
  12.10422696  19.46990488  31.31775346  50.37526829  81.02968365
 130.33795859 209.65136088] MeV
[1.72536301e+49 7.24984201e+49 3.63895427e+50 1.48246568e+51
 3.26752946e+51 2.92454494e+51 1.11961656e+51 1.47364614e+50
 1.52159412e+47 5.22545606e+42 9.13369778e+40 5.47718112e+39] erg / (MeV s)
[1.72536301e+49 7.24984201e+49 3.63895427e+50 1.48246568e+51
 3.26752946e+51 2.92454494e+51 1.11961656e+51 1.47364614e+50
 1.52159412e+47 5.22545606e+42 9.13369778e+40 5.47718112e+39] erg / (MeV s)
Max value:
7.525065547748337 MeV 3.267529456792245e+51 erg / (MeV s)

1
[ 1.07379938  1.57612195  2.31343064  3.3956518   4.98413523  7.31571003
 10.73799383 15.76121948 23.13430641 33.956518   49.84135225 73.15710034] MeV
[ 1.07379938  1.57612195  2.31343064  3.3956518   4.98413523  7.31571003
 10.73799383 15.76121948 23.13430641 33.956518   49.84135225 73.15710034] MeV
[5.92242985e+48 1.55280326e+49 4.50944197e+49 1.19821895e+50
 2.84295176e+50 5.81878888e+50 9.43454448e+50 1.10677342e+51
 7.89164432e+50 2.62666807e+50 4.37153391e+49 2.07419755e+47] erg / (MeV s)
[5.92242985e+48 1.55280326e+49 4.50944197e+49 1.19821895e+50
 2.84295176e+50 5.81878888e+50 9.43454448e+50 1.10677342e+51
 7.89164432e+50 2.62666807e+50 4.37153391e+49 2.07419755e+47] erg / (MeV s)
Max value:
15.761219478407002 MeV 1.1067734207970115e+51 erg / (MeV s)

2
[ 1.07379938  1.57612195  2.31343064  3.3956518   4.98413523  7.31571003
 10.73799383 15.76121948 23.13430641 33.956518   49.84135225 73.15710034] MeV
[ 1.07379938  1.57612195  2.31343064  3.3956518   4.98413523  7.31571003
 10.73799383 15.76121948 23.13430641 33.956518   49.84135225 73.15710034] MeV
[ 3.45572562e+48  1.29959365e+49  5.32565959e+49  2.04720220e+50
  7.01734472e+50  1.75827435e+51  2.77350985e+51  2.32469667e+51
  1.01137720e+51  2.01726857e+50  4.12624276e+48 -5.25592115e+43] erg / (MeV s)
[ 3.45572562e+48  1.29959365e+49  5.32565959e+49  2.04720220e+50
  7.01734472e+50  1.75827435e+51  2.77350985e+51  2.32469667e+51
  1.01137720e+51  2.01726857e+50  4.12624276e+48 -5.25592115e+43] erg / (MeV s)
Max value:
10.737993829320544 MeV 2.7735098481136573e+51 erg / (MeV s)

3
[ 1.07379938  1.57612195  2.31343064  3.3956518   4.98413523  7.31571003
 10.73799383 15.76121948 23.13430641 33.956518   49.84135225 73.15710034] MeV
[ 1.07379938  1.57612195  2.31343064  3.3956518   4.98413523  7.31571003
 10.73799383 15.76121948 23.13430641 33.956518   49.84135225 73.15710034] MeV
[5.92242985e+48 1.55280326e+49 4.50944197e+49 1.19821895e+50
 2.84295176e+50 5.81878888e+50 9.43454448e+50 1.10677342e+51
 7.89164432e+50 2.62666807e+50 4.37153391e+49 2.07419755e+47] erg / (MeV s)
[5.92242985e+48 1.55280326e+49 4.50944197e+49 1.19821895e+50
 2.84295176e+50 5.81878888e+50 9.43454448e+50 1.10677342e+51
 7.89164432e+50 2.62666807e+50 4.37153391e+49 2.07419755e+47] erg / (MeV s)
Max value:
15.761219478407002 MeV 1.1067734207970115e+51 erg / (MeV s)

[7]:

t = 600 * u.ms
E = np.arange(0,111) * u.MeV
theta = 23.55*u.degree
phi = 22.5*u.degree

spectra = model.get_initial_spectra(t, E, theta, phi)
cspectra = cmodel.get_initial_spectra(t, E, theta, phi)

fig, axes = plt.subplots(2,2, figsize=(11,8), sharex=True, sharey=True,
                         gridspec_kw = {'hspace':0.04, 'wspace':0.025})
axes = axes.flatten()

spectra = model.get_initial_spectra(t, E, theta, phi)

for ax, flavor in zip(axes, Flavor):
    ax.plot(E, spectra[flavor].to('1e52 erg/(s MeV)'), ls='-', lw=1, label='{}: model'.format(flavor.to_tex()))
    ax.plot(E, cspectra[flavor].to('1e52 erg/(s MeV)'), ls=':', lw=3, label='{}: cached model'.format(flavor.to_tex()))
    ax.legend(fontsize=14, ncol=1, loc='upper right')

    ax.set(xlim=(E[0].to_value('MeV'), E[-1].to_value('MeV')))
    ax.grid(ls=':')

axes[0].set(ylabel=r'$dL/dE$ [$10^{52}$ erg s$^{-1}$ MeV$^{-1}$]')
axes[2].set(xlabel=r'energy [MeV]')

fig.suptitle('Spectra at $t={{{:g}}}$'.format(t))
fig.subplots_adjust(top=0.925, bottom=0.1)
../../_images/nb_ccsn_Fornax_2019_11_0.png

Plot \(L_\nu(t,\theta,\varphi)\) at Several Times

Plot the angular dependence \(L(t)\) for several values of \(t\).

[8]:

if healpy_available:
    times = np.asarray([10., 100., 300., 600.]) * u.ms
    nt = len(times)

    fig, axes = plt.subplots(nt,3, figsize=(8*3, nt*4.5))

    L_nue = cmodel.luminosity[Flavor.NU_E].to('1e52*erg/s')
    L_nue_bar = cmodel.luminosity[Flavor.NU_E_BAR].to('1e52*erg/s')
    L_nux = cmodel.luminosity[Flavor.NU_X].to('1e52*erg/s')

    for i, t in enumerate(times):
        j = (np.abs(t - cmodel.time)).argmin()

        ax = axes[i,0]
        plt.axes(ax)
        hp.mollview(L_nue[j], title=r'$L_{{\nu_e}}(t={:.3f})$'.format(cmodel.time[j]),
                    unit=r'$10^{52}$ erg s$^{-1}$', cmap='viridis',
                    hold=True)

        ax = axes[i,1]
        plt.axes(ax)
        hp.mollview(L_nue_bar[j], title=r'$L_{{\bar{{\nu}}_e}}(t={:.3f})$'.format(cmodel.time[j]),
                    unit=r'$10^{52}$ erg s$^{-1}$', cmap='magma',
                    hold=True)

        ax = axes[i,2]
        plt.axes(ax)
        hp.mollview(L_nux[j], title=r'$L_{{\nu_X}}(t={:.3f})$'.format(cmodel.time[j]),
                    unit=r'$10^{52}$ erg s$^{-1}$', cmap='cividis',
                    hold=True)

Plot \(\Delta L_\nu(t,\theta,\varphi) / \langle L_\nu(t,\theta,\varphi)\rangle\) at Several Times

Same plot as above, but here plot the deviation from the average over the sky.

This demonstrates how the deviation from isotropy is very small initially and approaches 20% at later times.

[9]:

if healpy_available:
    times = np.asarray([10., 100., 300., 600.]) * u.ms
    nt = len(times)

    fig, axes = plt.subplots(nt,3, figsize=(8*3, nt*4.5))

    L_nue = cmodel.luminosity[Flavor.NU_E].to('1e52*erg/s')
    L_nue_bar = cmodel.luminosity[Flavor.NU_E_BAR].to('1e52*erg/s')
    L_nux = cmodel.luminosity[Flavor.NU_X].to('1e52*erg/s')

    vmin, vmax = -0.2, 0.2

    for i, t in enumerate(times):
        j = (np.abs(t - cmodel.time)).argmin()

        ax = axes[i,0]
        plt.axes(ax)
        Lavg = np.average(L_nue[j])
        dL_on_L = (L_nue[j] - Lavg) / Lavg
        hp.mollview(dL_on_L, title=r'$\nu_e$: $t=${:.3f}'.format(cmodel.time[j]),
                    unit=r'$(L_\nu - \langle L_\nu\rangle) / \langle L_\nu\rangle$', cmap='viridis',
                    min=vmin, max=vmax,
                    hold=True)

        ax = axes[i,1]
        plt.axes(ax)
        Lavg = np.average(L_nue_bar[j])
        dL_on_L = (L_nue_bar[j] - Lavg) / Lavg
        hp.mollview(dL_on_L, title=r'$\bar{{\nu}}_e$: $t=${:.3f}'.format(cmodel.time[j]),
                    unit=r'$\Delta L_{\bar{\nu}_e}/\langle L_{\bar{\nu}_e}\rangle$', cmap='magma',
                    min=vmin, max=vmax,
                    hold=True)

        ax = axes[i,2]
        plt.axes(ax)
        Lavg = np.average(L_nux[j])
        dL_on_L = (L_nux[j] - Lavg) / Lavg
        hp.mollview(dL_on_L, title=r'$\nu_X$: $t=${:.3f}'.format(cmodel.time[j]),
                    unit=r'$(L_\nu - \langle L_\nu\rangle) / \langle L_\nu\rangle$', cmap='cividis',
                    min=vmin, max=vmax,
                    hold=True)

Superimpose \(L_\nu(t,\theta,\varphi)\) at All Locations on the Sphere

Plot \(L(t)\) at all locations as well as the average, and then the deviation from the average.

Electron Neutrinos

[10]:

if healpy_available:
    Lavg = np.average(L_nue, axis=1)
    dL_over_L = (L_nue - Lavg[:, np.newaxis]) / Lavg[:, np.newaxis]

    fig, axes = plt.subplots(2, 1, figsize=(10, 8),
                             gridspec_kw={'height_ratios': [3, 1], 'hspace': 0},
                             sharex=True, tight_layout=True)

    ax = axes[0]
    ax.plot(cmodel.time, L_nue, color='gray', alpha=0.05)
    ax.plot(cmodel.time, Lavg, color='r', ls='--')
    ax.set(xlim=cmodel.time[0::len(cmodel.time)-1].value,
           ylabel=r'$L_{\nu_e}(t,\theta,\varphi)$ [$10^{52}$ erg s$^{-1}$]',
           ylim=(0, 1.1*np.max(L_nue).value),
           title=r'$L_{\nu_e}(t,\theta,\varphi)$')
    ax.grid(ls=':')

    ax = axes[1]
    ax.plot(cmodel.time, dL_over_L, color='gray', alpha=0.1)
    ax.plot(cmodel.time, np.zeros(cmodel.time.shape), ls='--', color='r')
    ax.set(xlabel='time [s]',
           ylabel=r'$(L_\nu - \bar{L}_\nu)/\bar{L}_\nu$',
           ylim=(-0.3, 0.3))
    ax.grid(ls=':')

Electron Antineutrinos

[11]:

if healpy_available:
    Lavg = np.average(L_nue_bar, axis=1)
    dL_over_L = (L_nue_bar - Lavg[:, np.newaxis]) / Lavg[:, np.newaxis]

    fig, axes = plt.subplots(2, 1, figsize=(10, 8),
                             gridspec_kw={'height_ratios': [3, 1], 'hspace': 0},
                             sharex=True, tight_layout=True)

    ax = axes[0]
    ax.plot(cmodel.time, L_nue_bar, color='gray', alpha=0.05)
    ax.plot(cmodel.time, Lavg, color='r', ls='--')
    ax.set(xlim=cmodel.time[0::len(cmodel.time)-1].value,
           ylabel=r'$L_{\bar{\nu}_e}(t,\theta,\varphi)$ [$10^{52}$ erg s$^{-1}$]',
           ylim=(0, 1.1*np.max(L_nue_bar).value),
           title=r'$L_{\bar{\nu}_e}(t,\theta,\varphi)$')
    ax.grid(ls=':')

    ax = axes[1]
    ax.plot(cmodel.time, dL_over_L, color='gray', alpha=0.1)
    ax.plot(cmodel.time, np.zeros(cmodel.time.shape), ls='--', color='r')
    ax.set(xlabel='time [s]',
           ylabel=r'$(L_\bar{\nu} - \bar{L}_\bar{\nu})/\bar{L}_\bar{\nu}$',
           ylim=(-0.3, 0.3))
    ax.grid(ls=':')

Flavor X (Stand-in for Mu and Tau Neutrinos)

Note that in this model, the mu and tau antineutrino flux is identical to the neutrino flux so we won’t bother plotting it here.

[12]:

if healpy_available:
    Lavg = np.average(L_nux, axis=1)
    dL_over_L = (L_nux - Lavg[:,np.newaxis]) / Lavg[:,np.newaxis]

    fig, axes = plt.subplots(2,1, figsize=(10,8),
                        gridspec_kw = {'height_ratios':[3,1], 'hspace':0},
                        sharex=True, tight_layout=True)

    ax = axes[0]
    ax.plot(cmodel.time, L_nux, color='gray', alpha=0.05)
    ax.plot(cmodel.time, Lavg, color='r', ls='--')
    ax.set(xlim=cmodel.time[0::len(cmodel.time)-1].value,
            ylabel=r'$L_{\nu_X}(t,\theta,\varphi)$ [$10^{52}$ erg s$^{-1}$]',
            ylim=(0, 1.1*np.max(L_nux).value),
            title=r'$L_{\nu_X}(t,\theta,\varphi)$')
    ax.grid(ls=':')

    ax = axes[1]
    ax.plot(cmodel.time, dL_over_L, color='gray', alpha=0.1)
    ax.plot(cmodel.time, np.zeros(cmodel.time.shape), ls='--', color='r')
    ax.set(xlabel='time [s]',
            ylabel=r'$(L_\nu - \bar{L}_\nu)/\bar{L}_\nu$',
            ylim=(-0.3, 0.3))
    ax.grid(ls=':');