{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Flavor Transformations\n", "\n", "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 " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from astropy import units as u \n", "\n", "from snewpy.neutrino import Flavor, MassHierarchy, MixingParameters\n", "from snewpy.models.ccsn import Nakazato_2013\n", "from snewpy.flavor_transformation import AdiabaticMSW, NonAdiabaticMSWH, \\\n", " TwoFlavorDecoherence, ThreeFlavorDecoherence, \\\n", " NeutrinoDecay, AdiabaticMSWes, NonAdiabaticMSWes\n", "\n", "mpl.rc('font', size=18)\n", "%matplotlib inline\n", "\n", "# `np.trapz` was renamed in numpy 2.0 and now raises a DeprecationWarning\n", "if np.__version__ < \"2.0\":\n", " np.trapezoid = np.trapz" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model = Nakazato_2013(progenitor_mass=20*u.solMass, revival_time=100*u.ms, metallicity=0.004, eos='shen')\n", "model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Transformed and untransformed energy integrated flux at Earth\n", "\n", "Compute and plot the flux at Earth for a SN at 10 kpc with no flavor transformation, and with the chosen flavor transformation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_total_flux(model, xform_nmo, xform_imo):\n", " \"\"\"Plot initial and oscillated neutrino luminosities.\n", " \n", " Parameters\n", " ----------\n", " model : SupernovaModel\n", " An input model from a CCSN simulation.\n", " flav_xform : FlavorTransformation\n", " A FlavorTransformation subclass; used to create an instance.\n", " \"\"\"\n", " \n", " energies = np.linspace(0,60,121) * u.MeV \n", " d = (10*u.kpc).to('cm').value # distance to SN\n", " \n", " times = model.get_time()\n", " burst_epoch = times <= 0.1*u.s\n", " accretion_epoch = (times > 0.1*u.s) & (times <= 0.5*u.s)\n", " cooling_epoch = (times > 0.5*u.s) & (times <= 10*u.s)\n", " \n", " ilum = {}\n", " olum_nmo = {}\n", " olum_imo = {}\n", " \n", " for flavor in Flavor:\n", " ilum[flavor] = np.zeros(len(times))\n", " olum_nmo[flavor] = np.zeros(len(times))\n", " olum_imo[flavor] = np.zeros(len(times))\n", "\n", " # Compute the transformed and untransformed flux at each time.\n", " for i, t in enumerate(times):\n", " ispec = model.get_initial_spectra(t, energies)\n", " ospec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)\n", " ospec_imo = model.get_transformed_spectra(t, energies, xform_imo)\n", "\n", " for flavor in Flavor:\n", " for j, E in enumerate(energies):\n", " ispec[flavor][j] /= (4.*np.pi*d**2)\n", " ospec_nmo[flavor][j] /= (4.*np.pi*d**2)\n", " ospec_imo[flavor][j] /= (4.*np.pi*d**2)\n", " \n", " for flavor in Flavor:\n", " ilum[flavor][i] = np.trapezoid(ispec[flavor].to('1/(erg*s)'), energies.to('erg')).value \n", " olum_nmo[flavor][i] = np.trapezoid(ospec_nmo[flavor].to('1/(erg*s)'), energies.to('erg')).value\n", " olum_imo[flavor][i] = np.trapezoid(ospec_imo[flavor].to('1/(erg*s)'), energies.to('erg')).value \n", " \n", " # make the figures \n", " fig, axes = plt.subplots(3,3, figsize=(20,12), tight_layout=True)\n", " \n", " smax = [0.,0.,0.]\n", " titles = ['Untransformed', 'Transformed (NMO)', 'Transformed (IMO)']\n", " for i, spec in enumerate([ilum, olum_nmo, olum_imo]):\n", " for j, phase in enumerate([burst_epoch, accretion_epoch, cooling_epoch]):\n", " ax = axes[i,j]\n", " timeunits = 'ms' if j==0 else 's'\n", " \n", " for flavor in Flavor:\n", " if i == 0:\n", " smax[j] = np.maximum(smax[j], 1.1*np.max(spec[flavor][phase]))\n", " \n", " ax.plot(times[phase].to(timeunits),\n", " spec[flavor][phase], label=flavor.to_tex(), lw=3,\n", " color='C0' if flavor.is_electron else 'C1',\n", " ls='-' if flavor.is_neutrino else ':')\n", " \n", " ax.set(xlim=(times[phase][0].to(timeunits).value, times[phase][-1].to(timeunits).value),\n", " ylim=(0, smax[j]))\n", " \n", " if j==0:\n", " ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')\n", " ax.legend(loc='upper right', ncol=1, fontsize=18)\n", " if j==1:\n", " ax.set(title=titles[i])\n", " if i < 2:\n", " ax.set(xticklabels=[])\n", " else:\n", " ax.set(xlabel='time [{}]'.format(timeunits))\n", " \n", " ax.grid(ls=':')\n", "\n", " \n", " return fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Untransformed and Transformed Spectra at Earth\n", "\n", "Compute and plot the spectra at Earth for a SN at 10 kpc with no flavor transformation, and with the chosen flavor transformation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_spectra(model, xform_nmo, xform_imo, t):\n", " \"\"\"Plot initial and oscillated neutrino luminosities.\n", " \n", " Parameters\n", " ----------\n", " model : SupernovaModel\n", " An input model from a CCSN simulation.\n", " flav_xform : FlavorTransformation\n", " A FlavorTransformation subclass; used to create an instance.\n", " t : astropy.Quantity\n", " Time to compute the spectrum.\n", " \"\"\"\n", "\n", " energies = np.linspace(0,60,121) * u.MeV \n", " d = (10*u.kpc).to('cm').value # distance to SN\n", "\n", " #get the spectra\n", " ispec = model.get_initial_spectra(t, energies) \n", " ospec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)\n", " ospec_imo = model.get_transformed_spectra(t, energies, xform_imo)\n", "\n", " for flavor in Flavor:\n", " for j, E in enumerate(energies):\n", " ispec[flavor][j] /= (4.*np.pi*d**2)\n", " ospec_nmo[flavor][j] /= (4.*np.pi*d**2)\n", " ospec_imo[flavor][j] /= (4.*np.pi*d**2)\n", " \n", " fig, axes = plt.subplots(2,2, figsize=(18,15), sharex=True, sharey=True, tight_layout=True)\n", "\n", " for i, spec in enumerate([ispec, ospec_nmo]):\n", " axes[0][0].plot(energies, spec[Flavor.NU_E]/1e16, \n", " label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),\n", " color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", " axes[0][0].plot(energies, spec[Flavor.NU_X]/1e16, \n", " label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),\n", " color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", "\n", " axes[0][0].set(title='Neutrinos in the NMO: $t = ${:.1f}'.format(t))\n", " axes[0][0].grid()\n", " axes[0][0].legend(loc='upper right', ncol=2, fontsize=16)\n", " \n", " axes[0][1].plot(energies, spec[Flavor.NU_E_BAR]/1e16, \n", " label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),\n", " color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", " axes[0][1].plot(energies, spec[Flavor.NU_X_BAR]/1e16, \n", " label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),\n", " color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", "\n", " axes[0][1].set(title='Antineutrinos in the NMO: $t = ${:.1f}'.format(t)) \n", " axes[0][1].grid()\n", " axes[0][1].legend(loc='upper right', ncol=2, fontsize=16) \n", " \n", " for i, spec in enumerate([ispec, ospec_imo]):\n", " axes[1][0].plot(energies, spec[Flavor.NU_E]/1e16, \n", " label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),\n", " color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", " axes[1][0].plot(energies, spec[Flavor.NU_X]/1e16, \n", " label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),\n", " color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", "\n", " axes[1][0].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Neutrinos in the IMO: $t = ${:.1f}'.format(t))\n", " axes[1][0].grid()\n", " axes[1][0].legend(loc='upper right', ncol=2, fontsize=16)\n", " \n", " axes[1][1].plot(energies, spec[Flavor.NU_E_BAR]/1e16, \n", " label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),\n", " color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", " axes[1][1].plot(energies, spec[Flavor.NU_X_BAR]/1e16, \n", " label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),\n", " color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", "\n", " axes[1][1].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Antineutrinos in the IMO: $t = ${:.1f}'.format(t))\n", " axes[1][1].grid()\n", " axes[1][1].legend(loc='upper right', ncol=2, fontsize=18) \n", "\n", " ax = axes[0][0]\n", " ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')\n", " ax = axes[1][0]\n", " ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')\n", " \n", " return fig" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plot_total_flux(model, AdiabaticMSW(), AdiabaticMSW(mh=MassHierarchy.INVERTED) )\n", "fig.show()\n", "# fig.savefig('flux_adiabaticmsw.pdf')\n", "fig = plot_spectra(model, AdiabaticMSW(), AdiabaticMSW(mh=MassHierarchy.INVERTED), 100*u.ms)\n", "fig.show()\n", "# fig.savefig('spectra_adiabaticmsw.pdf')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mp_nmo = MixingParameters()\n", "angles_nmo = mp_nmo.get_mixing_angles() + (10*u.deg,)\n", "xf_nmo = AdiabaticMSWes(angles_nmo)\n", "\n", "mp_imo = MixingParameters(mh=MassHierarchy.INVERTED)\n", "angles_imo = mp_imo.get_mixing_angles() + (10*u.deg,)\n", "xf_imo = AdiabaticMSWes(angles_imo)\n", "\n", "fig = plot_total_flux(model, xf_nmo, xf_imo)\n", "fig.show()\n", "# fig.savefig('flux_adiabaticmswes.pdf')\n", "fig = plot_spectra(model, xf_nmo, xf_imo, 100*u.ms)\n", "fig.show()\n", "# fig.savefig('spectra_adiabaticmswes.pdf')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.9.5 ('snews')", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" }, "vscode": { "interpreter": { "hash": "e2528887d751495e023d57d695389d9a04f4c4d2e5866aaf6dc03a1ed45c573e" } } }, "nbformat": 4, "nbformat_minor": 2 }