Kuroda 2020 Models

Models from Takami Kuroda. The paper describing the simulations is “Impact of magnetic field on neutrino-matter interactions in core-collapse supernova” arXiv:2009.07733. Included with SNEWPY with permission of the authors.

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

Kuroda_2020.param

{'rotational_velocity': <Quantity [0., 1.] rad / s>,
 'magnetic_field_exponent': [0, 12, 13],
 'eos': ['LS220'],
 'progenitor_mass': [<Quantity 20. solMass>]}

Quite a lot of choice there! However, for this 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 tuple of all combinations:
Kuroda_2020.get_param_combinations()

[{'rotational_velocity': <Quantity 0. rad / s>,
  'magnetic_field_exponent': 0,
  'eos': 'LS220',
  'progenitor_mass': <Quantity 20. solMass>},
 {'rotational_velocity': <Quantity 1. rad / s>,
  'magnetic_field_exponent': 12,
  'eos': 'LS220',
  'progenitor_mass': <Quantity 20. solMass>},
 {'rotational_velocity': <Quantity 1. rad / s>,
  'magnetic_field_exponent': 13,
  'eos': 'LS220',
  'progenitor_mass': <Quantity 20. solMass>}]

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

model = Kuroda_2020(**Kuroda_2020.get_param_combinations()[1])
# This is equivalent to:
# model = Kuroda_2020(rotational_velocity=1*u.rad/u.s, magnetic_field_exponent=12)
model

/tmp/ipykernel_1484/2294202880.py:1: UserWarning: Argument `eos` is deprecated.
  model = Kuroda_2020(**Kuroda_2020.get_param_combinations()[1])

Kuroda_2020 Model: LnuR10B12.dat

Parameter	Value
EOS	LS220
Progenitor mass	\(20\) \(\mathrm{M_{\odot}}\)
Rotational velocity	\(1\) \(\mathrm{\frac{rad}{s}}\)
B_0 Exponent	12

Finally, let’s plot the luminosity of different neutrino flavors for this model. (Note that the Kuroda_2020 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.0, 0.35),
       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 = 50*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();