Mori 2023 2D Models with Axionlike Production

Neutrino spectra from a set of 2D simulations with axionlike particle production. The models are described in the paper Multi-messenger signals of heavy axionlike particles in core-collapse supernovae: two-dimensional simulations by K. Mori, T. Takiwaki, K. Kotake and S. Horiuchi, Phys. Rev. D 108:063027, 2023.

The following models are supported:

Model	Axion mass [MeV]	Coupling g10	tpb,2000 [ms]	Ediag [1e51 erg]	MPNS [Msun]
Standard	−	0	390	0.40	1.78
(100, 2)	100	2	385	0.37	1.77
(100, 4)	100	4	362	0.34	1.76
(100, 10)	100	10	395	0.36	1.77
(100, 12)	100	12	357	0.43	1.77
(100, 14)	100	14	360	0.44	1.77
(100, 16)	100	16	367	0.51	1.77
(100, 20)	100	20	330	1.10	1.74
(200, 2)	200	2	374	0.45	1.77
(200, 4)	200	4	376	0.45	1.76
(200, 6)	200	6	333	0.54	1.75
(200, 8)	200	8	323	0.94	1.74
(200, 10)	200	10	319	1.61	1.73
(200, 20)	200	20	248	3.87	1.62

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

from astropy import units as u
from glob import glob

from scipy.interpolate import PchipInterpolator

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

mpl.rc('font', size=16)

Initialize the 2D models

Use the param property of the model class to see the available parameters. Models are initialized using the axion_mass and axion_coupling parameters.

Mori_2023.param

{'axion_mass': <Quantity [  0., 100., 200.] MeV>,
 'axion_coupling': <Quantity [ 0.,  2.,  4.,  6.,  8., 10., 12., 14., 16., 20.] 1e-10 / GeV>,
 'progenitor_mass': <Quantity [20.] solMass>}

The model with axion_mass=0 and axion_coupling=0 is a standard simulation with no ALP production.

model_std = Mori_2023(axion_mass=0, axion_coupling=0)
model_std.metadata

{'Progenitor mass': <Quantity 20. solMass>,
 'Axion mass': 0,
 'Axion coupling': <Quantity 0. 1e-10 / GeV>,
 'PNS mass': <Quantity 1.78 solMass>}

# Initialize a handful of axion models.
models = {}
for (am, ac) in ((100*u.MeV, 2e-10/u.GeV), (200*u.MeV,2e-10/u.GeV), (100*u.MeV,10e-10/u.GeV), (100*u.MeV,20e-10/u.GeV), (200*u.MeV,10e-10/u.GeV), (200*u.MeV,20e-10/u.GeV)):
    models[(am,ac)] = Mori_2023(axion_mass=am, axion_coupling=ac)

models

{(<Quantity 100. MeV>,
  <Quantity 2.e-10 1 / GeV>): Mori_2023 Model: t-prof_100_2.dat
 Progenitor mass  : 20.0 solMass
 Axion mass       : 100.0 MeV
 Axion coupling   : 2.0 1e-10 / GeV
 PNS mass         : 1.77 solMass,
 (<Quantity 200. MeV>,
  <Quantity 2.e-10 1 / GeV>): Mori_2023 Model: t-prof_200_2.dat
 Progenitor mass  : 20.0 solMass
 Axion mass       : 200.0 MeV
 Axion coupling   : 2.0 1e-10 / GeV
 PNS mass         : 1.77 solMass,
 (<Quantity 100. MeV>,
  <Quantity 1.e-09 1 / GeV>): Mori_2023 Model: t-prof_100_10.dat
 Progenitor mass  : 20.0 solMass
 Axion mass       : 100.0 MeV
 Axion coupling   : 10.0 1e-10 / GeV
 PNS mass         : 1.77 solMass,
 (<Quantity 100. MeV>,
  <Quantity 2.e-09 1 / GeV>): Mori_2023 Model: t-prof_100_20.dat
 Progenitor mass  : 20.0 solMass
 Axion mass       : 100.0 MeV
 Axion coupling   : 20.0 1e-10 / GeV
 PNS mass         : 1.74 solMass,
 (<Quantity 200. MeV>,
  <Quantity 1.e-09 1 / GeV>): Mori_2023 Model: t-prof_200_10.dat
 Progenitor mass  : 20.0 solMass
 Axion mass       : 200.0 MeV
 Axion coupling   : 10.0 1e-10 / GeV
 PNS mass         : 1.73 solMass,
 (<Quantity 200. MeV>,
  <Quantity 2.e-09 1 / GeV>): Mori_2023 Model: t-prof_200_20.dat
 Progenitor mass  : 20.0 solMass
 Axion mass       : 200.0 MeV
 Axion coupling   : 20.0 1e-10 / GeV
 PNS mass         : 1.62 solMass}

Plot Model Luminosities

Compare axion model luminosity to the standard 2D simulation.

Higher mass models with stronger coupling constants should produce a decrease in neutrino luminosity at all flavors relative to the reference simulation.

for (m,c), model in models.items():

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

    Lmin,  Lmax  = 1e99, -1e99
    dLmin, dLmax = 1e99, -1e99

    for j, (flavor) in enumerate([Flavor.NU_E, Flavor.NU_E_BAR, Flavor.NU_X, Flavor.NU_X_BAR]):
        ax = axes[0][j]

        ax.plot(model_std.time, model_std.luminosity[flavor]/1e51, 'k', label=r'$20M_\odot$ reference')  # Report luminosity in [foe/s]
        Lmin = np.minimum(Lmin, np.min(model_std.luminosity[flavor].to_value('1e51 erg/s')))
        Lmax = np.maximum(Lmax, np.max(model_std.luminosity[flavor].to_value('1e51 erg/s')))

        modlabel = rf"{flavor.to_tex()}: $m_a=${m.to_string(format='latex_inline')}\n     $g_{{a\gamma}}=${c.to_string(format='latex_inline')}"
        ax.plot(model.time, model.luminosity[flavor]/1e51,  # Report luminosity in units foe/s
                label=modlabel,
                color='C0' if flavor.is_electron else 'C1',
                ls='-' if flavor.is_neutrino else ':',
                lw=2)
        if j==0:
            ax.set(ylabel=r'luminosity [$10^{51}$ erg s$^{-1}$]')

        ax.legend(fontsize=12)
        ax.set(xlim=(model_std.time[0].to_value('s'), model_std.time[-1].to_value('s')))

        ax = axes[1][j]
        tmin = np.maximum(model.time[0], model_std.time[0]).to_value('s')
        tmax = np.minimum(model.time[-1], model_std.time[-1]).to_value('s')
        times = np.arange(tmin, tmax, 0.001)*u.s

        Lstd = PchipInterpolator(model_std.time, model_std.luminosity[flavor].to_value('1e51 erg/s'))
        Lstd_t = Lstd(times)
        select = Lstd_t != 0

        Lmod = PchipInterpolator(model.time, model.luminosity[flavor].to_value('1e51 erg/s'))
        Lmod_t = Lmod(times)
        dL = (Lmod_t[select] - Lstd_t[select]) / Lstd_t[select]

        dLmin = np.minimum(dLmin, np.min(dL))
        dLmax = np.maximum(dLmax, np.max(dL))

        ax.plot(times[select], dL)
        if j==0:
            ax.set(xlabel='time [s]',
                   ylabel=r'$\Delta L_\nu/L_\nu$')

    for j in range(4):
        axes[0][j].set(ylim=(Lmin, 1.1*Lmax))
        axes[1][j].set(ylim=(dLmin, dLmax))

    fig.suptitle(rf"Axionlike model: $m_a=${m.to_string(format='latex_inline')}, $g_{{a\gamma}}=${c.to_string(format='latex_inline')}")